Army  Research  Laboratory 


Theoretical  Studies  of 
Solid  Nitromethane 


by  Dan  C.  Sorescu,  Betsy  M.  Rice, 
and  Donald  L.  Thompson 


ARL-RP-14 


February  2001 


A  reprint  from  the  Journal  of  Physical  Chemistry  B ,  vol.  104,  no.  35,  pp.  8406-8419. 


20010316  069 


Approved  for  public  release;  distribution  is  unlimited. 


The  findings  in  this  report  are  not  to  be  construed  as  an  official 
Department  of  the  Army  position  unless  so  designated  by  other 
authorized  documents. 

Citation  of  manufacturer’s  or  trade  names  does  not  constitute  an 
official  endorsement  or  approval  of  the  use  thereof. 

Destroy  this  report  when  it  is  no  longer  needed.  Do  not  return 
it  to  the  originator. 


Army  Research  Laboratory 

Aberdeen  Proving  Ground,  MD  21005-5066 


ARL-RP-14  _ February  2001 


Theoretical  Studies  of  Solid  Nitromethane 


Dan  C.  Sorescu 

Oklahoma  State  University 

Betsy  M.  Rice 

Weapons  and  Materials  Research  Directorate,  ARL 

Donald  L.  Thompson 

Oklahoma  State  University 


A  reprint  from  the  Journal  of  Physical  Chemistry  B,  vol.  104,  no.  35,  pp.  8406-8419. 


Approved  for  public  release;  distribution  is  unlimited. 


Abstract 


A  classical  potential  to  simulate  the  dynamics  of  a  nitromethane  crystal  as  a  function 
of  temperature  and  pressure  is  described.  The  intramolecular  part  of  the  potential  was 
taken  as  superposition  of  bond  stretching,  bond  bending,  and  torsional  angles  terms. 
These  terms  were  parametrized  on  the  basis  of  the  geometric  and  spectroscopic 
(vibrational  frequencies  and  eigenvectors)  data  obtained  using  ab  initio  molecular  orbital 
calculations  performed  at  the  B3LYP/6-31G  level  on  an  isolated  molecule.  The 
intermolecular  potential  used  is  of  the  Buckingham  6-exp  form  plus  charge-charge 
Coulombic  interactions  and  has  been  previously  developed  by  us  (Sorescu,  D.  C.;  Rice, 
B.  M.;  Thompson,  D.  L.  J.  Phys.  Chem.  1997,  B101,  798)  to  simulate  crystals  containing 
nitramine  molecules  and  several  other  classes  of  nitro  compounds.  The  analyses 
performed  using  constant  pressure  and  temperature  molecular  dynamics  simulations  and 
molecular  packing  calculations  indicate  that  the  proposed  potential  model  is  able  to 
reproduce  accurately  the  changes  of  the  structural  crystallographic  parameters  as 
functions  of  temperature  or  pressure  for  the  entire  range  of  values  investigated.  In 
addition,  the  calculated  bulk  modulus  of  nitromethane  was  found  in  excellent  agreement 
with  the  corresponding  experimental  results.  Moreover,  it  was  determined  that  the 
present  potential  predicts  correctly  an  experimentally  observed  45°  change  in  methyl 
group  orientation  in  the  high-pressure  regime  relative  to  the  low-temperature 
configuration.  The  analysis  of  the  linear  expansion  coefficients  and  linear  compression 
data  indicate  anisotropic  behavior  for  the  unit  cell  edges. 
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A  classical  potential  to  simulate  the  dynamics  of  a  nitromethane  crystal  as  a  function  of  temperature  and 
pressure  is  described.  The  intramolecular  part  of  the  potential  was  taken  as  superposition  of  bond  stretching, 
bond  bending,  and  torsional  angles  terms.  These  terms  were  parametrized  on  the  basis  of  the  geometric  and 
spectroscopic  (vibrational  frequencies  and  eigenvectors)  data  obtained  using  ab  initio  molecular  orbital 
calculations  performed  at  the  B3LYP/6-31G*  level  on  an  isolated  molecule.  The  intermolecular  potential 
used  is  of  the  Buckingham  6-exp  form  plus  charge— charge  Coulombic  interactions  and  has  been  previously 
developed  by  us  (Sorescu,  D.  C.;  Rice,  B.  M.;  Thompson,  D.  L.  J.  Phys.  Chem.  1997,  B101,  798)  to  simulate 
crystals  containing  nitramine  molecules  and  several  other  classes  of  nitro  compounds.  The  analyses  performed 
using  constant  pressure  and  temperature  molecular  dynamics  simulations  and  molecular  packing  calculations 
indicate  that  the  proposed  potential  model  is  able  to  reproduce  accurately  the  changes  of  the  structural 
crystallographic  parameters  as  functions  of  temperature  or  pressure  for  the  entire  range  of  values  investigated. 
In  addition,  the  calculated  bulk  modulus  of  nitromethane  was  found  in  excellent  agreement  with  the 
corresponding  experimental  results.  Moreover,  it  was  determined  that  the  present  potential  predicts  correctly 
an  experimentally  observed  45°  change  in  methyl  group  orientation  in  the  high-pressure  regime  relative  to 
the  low-temperature  configuration.  The  analysis  of  the  linear  expansion  coefficients  and  linear  compression 
data  indicate  anisotropic  behavior  for  the  unit  cell  edges. 


I.  Introduction 

This  is  the  seventh  in  a  series  of  papers  describing  our 
development  and  assessment  of  interaction  potentials  to  be  used 
in  the  study  of  dynamic  processes  in  energetic  materials.  Our 
original  intent  was  to  develop  a  model  to  study  nonreactive 
processes  in  the  nitramine  explosive  RDX  (1,3,5-hexahydro- 
1, 3,5-5- triazine).1  The  potential  that  was  developed  consisted 
of  atom— atom  (6-exp)  Buckingham  potential  terms  plus  elec¬ 
trostatic  interactions.  The  Coulombic  interactions  were  obtained 
through  fitting  atom-centered  partial  charges  to  a  quantum- 
mechanically-determined  electrostatic  potential  for  a  single  RDX 
molecule  whose  structure  corresponded  to  that  in  the  crystal  at 
ambient  conditions.  The  remaining  Buckingham  parameters 
were  adjusted  to  reproduce  the  experimental  structure  of  the 
RDX  crystal  at  ambient  conditions.  We  found  that  this  interac¬ 
tion  potential  could  also  describe  the  geometric  parameters  and 
lattice  energies  of  different  polymorphic  phases  of  two  other 
nitramine  crystals:  the  polycyclic  nitramine  2,4,6,8,10,12- 
hexanitrohexaazaisowurtzitane  (HNTW,  or  CL-20)2  and  the 
monocyclic  nitramine  octahydro-l,3,5,7-tetranitro-l,3,5,7-tet- 
raazacyclooctane  (HMX).3  Isothermal— isobaric  molecular  dy¬ 
namics  (NPT-MD)  simulations  for  these  crystals  predicted  cell 
parameters  within  a  few  percent  of  experimental  values,  and 
little  translational  or  rotational  disorder  of  the  molecules.  Further 
investigations  to  explore  the  limits  of  transferability  of  this 
interaction  potential  to  other  energetic  molecular  crystals  were 
undertaken  through  performing  molecular  packing  (MP) 
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calculations  for  30  nitramine  crystals.  These  included  several 
types  of  mono-  and  polycyclic  nitramines,  particularly  crystals 
of  importance  in  energetic  materials.4  For  most  of  the  crystals, 
the  predicted  structural  lattice  parameters  obtained  from  mo¬ 
lecular  packing  (MP)  calculations  differ  by  less  than  2%  from 
the  experimental  structures,  with  small  rotations  and  practically 
no  translations  of  the  molecules  in  the  asymmetric  unit  cell.  A 
further  assessment  on  the  limits  of  transferability  of  the 
interaction  potential  was  accomplished  through  molecular 
packing  calculations  of  51  crystals  containing  non- nitramine 
molecules  with  functional  groups  common  to  energetic  materi¬ 
als.5  MP  calculations  using  this  interaction  potential  reproduced 
the  crystal  structures  to  within  5%  of  experiment  for  these  51 
non-nitramine  systems,  including  the  explosives  pentaerythritol 
tetranitrate  (PETN),  nitromethane,  2,4,6-trinitrotoluene  (TNT), 
and  several  nitrocubane  derivatives. 

More  recently  we  have  analyzed  the  dynamics  of  the 
important  energetic  crystals  RDX,  HMX,  HNTW,  and  PETN 
under  hydrostatic  compression  conditions  using  NPT-MD 
simulations  and  this  intermolecular  potential.6  In  that  study  we 
found  that  the  predicted  lattice  parameters  for  the  RDX,  HMX, 
and  HNTW  crystals  are  in  good  agreement  with  experimental 
values  over  the  entire  range  of  pressures  investigated  experi¬ 
mentally.  For  the  PETN  crystal,  the  calculated  crystallographic 
parameters  were  in  acceptable  agreement  with  experimental  data 
for  pressures  up  to  5  GPa.  For  higher  pressures,  the  disagree¬ 
ments  of  predictions  and  experiment  were  attributed  to  the 
inadequacy  of  the  rigid-body  approximation  used  to  simulate 
floppy  molecules  such  as  PETN. 

The  validity  of  the  rigid  molecular  approximation  was 
assumed  in  our  previous  studies.1-6  As  indicated  by  our  results, 
this  model  has  been  very  successful  in  its  ability  to  describe 
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Solid  Nitromethane 


Figure  1.  Representation  of  the  nitromethane  crystal  unit  cell  with 
orthorhombic  space  group  P2i2i2i  and  Z  =  4  molecules  per  unit  cell. 
Atom  labels  are  consistent  with  the  indices  given  in  Table  1. 


the  equilibrium  structures  of  a  variety  of  organic  molecular 
crystals  under  ambient  conditions  and  with  moderate  increases 
in  pressure  and  temperature.  However,  the  physical  and  chemical 
processes  of  energetic  materials  that  are  of  most  interest  occur 
within  the  regime  of  high  pressures  and  temperatures,  a  regime 
in  which  conformational  molecular  changes  become  important. 
Consequently,  further  developments  of  the  interaction  potential 
are  necessary  to  describe  more  realistically  the  intramolecular 
motion,  molecular  deformations,  and  the  energy  flow  inside 
these  crystals.  For  this  purpose,  in  the  present  paper  we  eliminate 
the  previously  used  rigid-molecule  approximation1"6  and  extend 
the  current  intermolecular  potential  to  include  a  full  intra¬ 
molecular  potential  for  use  in  simulations  of  energetic  materials. 
Particularly,  in  this  work  we  consider  the  prototypical  explosive, 
nitromethane. 

Nitromethane  was  selected  for  our  first  attempt  at  developing 
a  fully  flexible  model  of  an  energetic  molecular  crystal  since 
numerous  experimental  investigations  of  its  properties  under  a 
wide  range  of  conditions  have  been  performed,  thus  providing 
significant  data  for  use  in  assessing  the  model  potential.7"14 
Single-crystal  X-ray  diffraction  and  neutron  powder  diffraction 
measurements  indicate  that  the  low-temperature  structure  of 
crystalline  nitromethane  belongs  to  the  orthorhombic  space 
group  P2i2i2i  with  Z  =  4  molecules  in  the  unit  cell  (see  Figure 
l).7  In  this  configuration  the  methyl  group  is  in  the  staggered 
symmetric  position  with  respect  to  the  C— NO2  plane.  At 
ambient  pressure,  the  crystal  symmetry  was  observed  to  remain 
unchanged  when  the  temperature  was  increased  from  4.2  to  228 
K.7’9  In  the  temperature  range  243—244  K  an  intermediate 
crystal  phase  (mesophase)  was  recently  observed,10  followed 
by  transition  to  the  liquid  state  at  244.73  K.11  No  crystal¬ 
lographic  details  have  been  provided  about  the  intermediate 
phase  at  243  K.10 

The  internal  rotation  of  the  methyl  group  in  the  crystal  has 
received  significant  attention,  since  experimental  evidence 
indicates  that  this  motion  is  almost  completely  governed  by 
intermolecular  interactions.  Analyses  of  the  neutron  powder 
diffraction  patterns7  for  deuterated  nitromethane  at  temperatures 
ranging  from  25  to  125  K  provided  temperature-dependent 
amplitudes  of  librations  of  the  methyl  group  about  the  C-N 
bond.  Quasielastic  neutron  scattering  spectra8  recorded  at 
temperatures  ranging  from  50  to  150  K  are  consistent  with 
rotations  of  the  methyl  groups  through  120°  jumps  rather  than 
the  60°  jumps  observed  in  the  gas  phase.  Furthermore,  it  was 
found8  that  in  the  crystalline  phase  the  activation  energy  for 
the  internal  rotation  of  the  methyl  group  is  about  234  cal/mol, 
significantly  higher  than  the  corresponding  gas-phase  barrier 
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of  6  cal/mol.15  Inelastic  neutron-scattering  measurements14  of 
tunneling  splittings  and  of  the  lowest  rotational  states  indicated 
that  the  methyl  group  activation  barrier  is  significantly  smaller 
than  the  energy  of  the  second  excited  rotational  state  by  ~  170 
cal/mol.  These  findings  have  been  confirmed  by  infrared16  and 
Raman10  measurements,  but  the  corresponding  activation  ener¬ 
gies  reported  in  these  studies  are  slightly  higher,  namely  401 
and  576  cal/mol,  respectively.  Several  mechanisms  have  been 
proposed  to  contribute  to  the  dynamics  of  methyl  groups 
including  the  enhanced  hopping  due  to  interactions  with 
phonons,  band  tunneling  through  thermally  broadened  excited 
states,  or  changes  of  the  barrier  with  temperature.14 

The  influence  of  pressure  on  the  internal  rotational  motion 
of  the  methyl  group  has  also  been  examined.  The  crystal 
structure  for  isotropically  compressed  nitromethane  has  been 
determined  through  X-ray  diffraction' at  pressures  ranging  from 
0.3  to  15.3  GPa,  at  room  temperature.9’12  It  was  found  that  the 
space  group  and  the  packing  arrangement  remain  unchanged 
relative  to  the  low-temperature  structure.7  However,  the  mea¬ 
surements  indicate  a  continuous  change  in  the  orientation  of 
the  methyl  group  with  the  increase  of  pressure.  For  pressures 
below  0.6  GPa  the  methyl  group  is  freely  rotating.  At  intermedi¬ 
ate  pressures  (between  0.6  and  3.5  GPa)  the  rotation  of  the 
methyl  group  is  hindered,  while  for  pressures  above  3.5  GPa 
the  orientation  of  the  methyl  group  becomes  fixed.  Moreover, 
it  was  determined  that  at  pressures  above  3.5  GPa,  the  methyl 
group  is  rotated  by  45°  relative  to  the  low-temperature  config¬ 
uration.7 

The  ensemble  of  the  above-presented  experimental  data 
provides  many  metrics  against  which  to  assess  the  performance 
of  our  interaction  potential  by  molecular  dynamics  simulations. 
In  this  work,  we  describe  such  an  assessment.  In  Section  n,  we 
provide  a  description  of  our  potential  model  and  our  choice  and 
parametrization  of  the  intramolecular  and  intermolecular  interac¬ 
tion  terms.  In  Section  HI,  we  provide  a  brief  description  of 
computational  details.  Section  IV  contains  the  results  of  our 
MP  and  NPT-MD  simulations  and  a  comparison  of  predictions 
of  crystal  structural  parameters  with  experimental  values  over 
a  wide  range  of  temperatures  and  pressures.  Summary  and 
conclusions  are  given  in  Section  V. 

II.  Potential  Energy  Functions 

We  have  assumed  that  the  potential  energy  for  a  system  of 
N  nitromethane  molecules  can  be  described  as  the  sum  of  inter- 
and  intramolecular  interaction  terms: 

N  I  j  N 

ytotal  __  j  ^intramolecular  _|_  _ 

i=l  \  2j=i 

The  intermolecular  potential  is  the  same  as  described  in  the 
earlier  studies,1"6  and  consists  of  the  superposition  of  a  pairwise 
sum  of  Buckingham  (6-exp)  (repulsion  and  dispersion)  and 
coulombic  (C)  potentials  of  the  form: 

Vofi(r)  =  A0fi  exp(-5^r)  -  C^r  (2) 


ij 


itermolecular  I 


(i) 


and 


vJ(r)  = 


4jt e0r 


(3) 


where  r  is  the  interatomic  distance  between  atoms  a  and  /?,  qa 
and  qp  are  the  electrostatic  charges  on  the  atoms,  and  <=o  is  the 
dielectric  permittivity  constant  of  free  space.  The  parameters 
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TABLE  1:  Force  Field  Parameters  for  Crystalline 
Nitromethane'1 


atom* 

charge/e 

c, 

-0.305342 

n2 

0.820603 

03 

—0.470445 

04 

-0.484277 

Hs 

0.143365 

H« 

0.155443 

H, 

0.140652 

Morse  potential  parameters 

bond 

(kJ/mol)  /S(A-') 

M) 

C-N 

N-O 

C-H 

251.0419414  2.005501 

390.3702976  2.459942 

426.7713101  1.892486 

1.499574 

1.226747 

1.090000 

bending  potential  parameters 

angle 

WkJ/mol  rad'2) 

0(deg) 

C-N-O 

O-N-O 

N-C-H 

H-C-H 

294.5211162 

657.8485114 

224.8140504 

149.9402928 

117.039960 

125.890000 

107.560708 

111.312289 

torsional  potential  parameters 

angle 

V*  (kJ/mol) 

<5(deg)  m 

HrC,  -N2-O3  (i  =  5,6,7)  0.27000  (0.49)c 

HrCj  -N2-O4  (1  =  5,6,7)  0.27000  (0.49)c 

N2-O4-O3-C1  240.37076 

-90.0  3.0 

90.0  3.0 

180.0  2.0 

a  The  electrostatic  charges  have  been  determined  by  the  CHELPG 
procedure  as  implemented  in  Gaussian  94 17  at  the  HF/6-31G**  level. 
*The  atom  designation  numbers  are  defined  in  Figure  1.  c  The  choice 
of  a  torsional  barrier  of  0.49  kJ/mol  is  also  discussed  in  text. 

Aofiy  and  for  different  types  of  atomic  pairs  have  been 
previously  published  and  have  been  used  in  the  present  study 
without  change.1 

The  set  of  partial  charges  used  in  these  calculations  were 
determined  through  fitting  these  to  the  quantum-mechanically- 
derived  electrostatic  interaction  potential  for  an  isolated  molecule 
whose  atoms  are  arranged  in  the  experimental  crystallographic 
arrangement.  These  calculations  have  been  done  using  the 
CHELPG  procedure  as  implemented  in  the  Gmissian  94 
package.17  We  have  previously  shown45  that  for  the  majority 
of  the  crystals  studied,  the  best  agreement  between  simulations 
and  experiment  occurs  when  the  set  of  partial  charges  is 
determined  using  methods  that  employ  electron  correlation 
effects  such  as  second-order  Moller-Plesset  (MP2)  perturbation 
theory.18  However,  this  was  not  true  for  nitromethane;  the  best 
agreement  between  the  measured  and  the  predicted  crystal¬ 
lographic  lattice  parameters  was  obtained  when  unsealed  partial 
charges  derived  from  the  Hartree— Fock  (HF)  wave  function 
were  used.  The  lack  of  improvement  in  the  predicted  lattice 
parameters  when  electron  correlation  effects  are  considered 
might  be  due  to  omission  of  other  important  interactions.  These 
could  include  the  omission  of  the  polarization  effects  of 
neighboring  molecules  in  the  crystal  when  evaluating  the 
electrostatic  charges  for  the  isolated  molecule,  or  the  assumption 
of  a  simple  charge -charge  interaction  model  when  a  multipolar 
interaction  description  would  be  more  accurate.  Nevertheless, 
our  choice  for  the  set  of  charges  determined  at  the  HF  level 
(see  Table  1),  together  with  the  rest  of  the  potential  parameters, 
give  a  very  good  representation  of  the  crystallographic  param¬ 
eters. 

In  a  previous  study  we  have  developed  a  complex  intramo¬ 
lecular  interaction  potential  to  describe  the  unimolecular 
decomposition  reactions  of  gas-phase  nitromethane.19  However, 


we  do  not  require  such  a  complex  function  in  order  to 
accomplish  the  goals  described  in  this  work.  Thus,  we  have 
assumed  a  simpler  form  of  the  intramolecular  interaction 
potential: 


yini 


tramolecular _ 


y  4-1/  4-  y 

stretch  v  bend  torsion 


(4) 


to  describe  the  bond  stretching,  angle  bending,  and  torsional 
motions  that  occur  within  an  isolated  molecule.  The  bond 
stretches  are  represented  by  Morse  potentials, 

6 

vstrctch  =  X  OJexp[-2A(ri  -  r,0)]  -  2  exp[-ft(r,  -  r,0)]} 

(5) 

where  r,  are  the  bond  distances,  Dc [  are  the  bond  dissociation 
energies,  /?,  are  the  curvature  parameters,  and  n°  are  the 
equilibrium  bond  lengths.  The  six  bond  stretches  that  are 
described  by  this  potential  correspond  to  the  six  covalent  bonds 
in  nitromethane. 

The  bending  potentials  are  represented  by  harmonic  functions 
of  the  form 


9  1 

V,bend  =  I;W-ei°)2  (6) 

where  k&  is  the  force  constant  and  0°  the  equilibrium  value  of 
the  angle.  The  nine  angles  used  to  describe  these  interactions 
correspond  to  the  HCH,  HCN,  CNO,  and  ONO  angles  in 
nitromethane. 

Cosine-type  torsional  potentials  have  been  used  to  describe 
the  relative  positions  of  the  C— NO2  atoms  and  the  orientations 
of  the  H  atoms  relative  to  the  C-N-0  planes.  The  potentials 
have  taken  the  form 


7 

torsion  =  X  +  C°S( ~  6 )]  (7) 

1=1 

where  V#  is  half  of  the  intramolecular  torsional  barrier,  is 
the  torsional  angle,  and  m  =  2  or  3.  Six  of  the  torsional  angles 
used  to  describe  these  interactions  correspond  to  HCNO  angles, 
the  seventh  corresponds  to  the  N-O-O-C  torsional  angle.  The 
complete  list  of  potential  parameters  from  eqs  5—7  is  given  in 
Table  1. 

The  bond  dissociation  energies  in  eq  5  have  been  taken  from 
our  previous  study  of  nitromethane  in  the  gas  phase,19  while 
the  remaining  set  of  geometrical  equilibrium  values  and  force 
constants  in  eqs  5—7  were  parametrized  based  on  quantum 
mechanical  information  generated  using  density  functional 
theory.  These  calculations  have  been  done  for  the  isolated 
molecule  using  Becke’s  three-parameter  hybrid  method,203  in 
combination  with  the  Lee,  Yang,  and  Pair  correlation  functional206 
(B3LYP)  and  the  basis  set  6-3 1G*21  using  the  Gaussian  94 
package  of  programs.17  Our  earlier  molecular  dynamics  studies22 
on  reactions  of  gas-phase  energetic  molecules  suggested  that  a 
potential  energy  function  fitted  to  only  the  structural  parameters 
and  normal-mode  frequencies  could  lead  to  a  model  that 
produced  anomalous  and  unexpected  results,  such  as  nonstatis- 
tical  behavior  in  the  unimolecular  decomposition  of  the  mol¬ 
ecule.  However,  when  the  function  was  fitted  to  the  ab  initio 
Cartesian  second  derivatives  of  the  energy,  the  potential  model 
produced  the  expected  statistical  behavior  in  the  unimolecular 
decomposition  of  the  molecule.  In  this  study,  both  the  geo¬ 
metrical  parameters  of  the  optimized  nitromethane  molecule  and 
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TABLE  2:  Comparison  between  the  Experimental,  the  Scaled  Harmonic  Vibrational  Frequencies  Calculated  at  B3LYP/ 
6-31G*LeveI,a  and  the  Predicted  Frequencies  by  Our  Force  Fields;  the  Corresponding  Projections  (Scalar  Product)  of  the 
Predicted  Eigenvectors  on  the  Ab  Initio  Eigenvectors  Are  Also  Indicated*  _ 


band 

type 

ref 

24a 

ref 

24b 

ref 

24c 

ref 

24d 

ref 

24e 

B3LYP/ 

6-31G*fl 

Mgas 

Pgas 

Ml 

PI 

M2 

P2 

Vt 

51 

50 

0.94 

117 

0.94 

157 

0.94 

Vs 

A" 

476 

477 

480 

475.2 

461 

459 

0.99 

461 

0.99 

461 

0.99 

V6 

4 

599 

605 

607 

602.5 

584 

596 

0.94 

601 

0.94 

605 

0.94 

vs 

A' 

647 

658 

656.0 

655 

657.4 

636 

602 

0.92 

605 

0.92 

605 

0.92 

V4 

A' 

921 

918 

917.1 

917 

917.9 

892 

828 

0.93 

821 

0.92 

821 

0.92 

V\\ 

ET 

1097 

1087 

1091.0 

1103 

1099.0 

1076 

1052 

0.99 

1053 

0.99 

1054 

0.92 

Vn 

ET 

1153 

1100 

1146.0 

1103 

1119.0 

1101 

1098 

0.99 

1101 

0.99 

1103 

0.99 

V3 

A ' 

1384 

1377 

1378.8 

1379 

1380.4 

1364 

1509 

0.95 

1508 

0.95 

1508 

0.99 

V2 

A' 

1397 

1397.0 

1402 

1397.4 

1388 

1372 

0.96 

1368 

0.96 

1368 

0.95 

Vio 

Ef 

1449 

1443 

1426 

1432 

1432 

1.00 

1434 

1.00 

1437 

0.99 

VlO 

ET 

1488 

1482 

1438.0 

1426 

1444 

1434 

0.99 

1436 

0.99 

1439 

0.99 

Vl 

A'4 

1582 

1586 

1583.3 

1561 

1583.8 

1615 

1647 

1.00 

1618 

1.00 

1618 

1.00 

V\ 

A' 

2965 

2972 

2964.3 

2967 

2973.9 

2981 

2961 

1.00 

2961 

0.98 

2961 

0.98 

v9 

ET 

2984.0 

3045 

3044.0 

3070 

3093 

0.98 

3093 

1.00 

3093 

0.98 

v9 

ET 

3048 

3065 

3080.4 

3065 

3080.0 

3100 

3095 

0.98 

3095 

1.00 

3095 

1.00 

« B3LYP/6-31G*  values  are  scaled  by  0.9613,  as  recommended  in  ref  23.  *The  columns  denoted  with  Mgas,  Ml,  and  M2  correspond  to  the 
reference  force  fields  for  gas  phase,  for  solid  phase  (with  torsional  HCNO  barrier  of  0.27  kJ/mol),  and  for  solid  phase  with  increased  torsional 
barrier  (torsional  HCNO  barrier  of  0.49  kJ/mol),  respectively.  The  corresponding  projections  of  the  predicted  eigenvectors  are  indicated  in  columns 
Pgas,  PI,  and  P2,  respectively.  The  frequency  units  are  cm"1. 


the  ab  initio  eigenvalues  (scaled  by  a  factor  of  0.96 1323)  and 
the  corresponding  eigenvectors  of  the  ab  initio  Hessian  matrix 
have  been  used  in  the  fitting  procedure. 

Special  attention  was  paid  in  this  work  to  the  choice  of  the 
torsional  barrier  V$>  corresponding  to  the  H-C— N— O  torsional 
angles.  As  indicated  in  previous  experimental  studies,8’101416 
an  exact  value  for  this  barrier  is  not  yet  well-known  and  several 
values  have  been  proposed.  The  major  difficulties  in  determining 
a  precise  barrier  for  the  internal  rotation  of  the  methyl  group 
were  due  to  the  existence  of  a  complex  mechanism  that 
incorporates  enhanced  hopping  due  to  interactions  with  phonons, 
tunneling  effects,  and  changes  of  the  barrier  with  temperature. 
However,  it  is  important  to  point  out  that  independent  of 
experimental  technique  used,  the  numerical  values  of  the 
activation  energy  for  the  methyl  rotation  are  quite  small81016 
and  represent  only  fractions  of  1  kcal/mol.  In  the  present  study 
we  have  investigated  several  values  for  V<x>  potential  parameters 
of  H-C— N— O  torsions  and  have  analyzed  their  role  on  the 
crystallographic  parameters.  The  vibrational  frequencies  pre¬ 
dicted  by  our  classical  potential  for  an  isolated  nitromethane 
molecule  are  given  in  Table  2  along  with  the  corresponding  ab 
initio  values  and  experimental  data.  In  this  table  we  also  indicate 
the  values  of  the  projections  of  the  eigenvectors  obtained  from 
the  normal-mode  analysis  using  the  proposed  potential  onto  their 
quantum  mechanical  counterpart.  A  projection  whose  magnitude 
is  one  would  indicate  perfect  agreement  between  the  model  and 
the  B3LYP/6-31G*  eigenvectors.  A  projection  that  has  a 
magnitude  near  zero  indicates  that  the  atomic  motion  of  the 
vibration  that  is  predicted  by  the  model  is  extremely  different 
from  that  predicted  by  the  quantum  mechanical  calculations. 
The  projections  of  the  normal-mode  eigenvectors  generated  by 
different  potential  models  (see  below)  onto  the  quantum 
mechanical  eigenvectors  are  also  indicated  in  Table  2. 

In  an  initial  attempt  we  described  the  internal  methyl  group 
rotational  motion  using  a  6-fold  torsional  potential  with  a  barrier 
of  6  cal/mol,  in  agreement  with  the  experimental  gas-phase 
barrier.15  The  vibrational  frequencies  corresponding  to  this 
model  are  indicated  in  Table  2  where  they  are  denoted  as  Mgas. 
As  can  be  seen  there  is  an  overall  very  good  agreement  between 
the  entire  set  of  predicted  vibrational  frequencies  and  the 
corresponding  ab  initio  results.  In  particular,  the  frequency  of 
the  torsional  mode  vT  is  practically  identical  to  the  corresponding 


quantum  mechanical  value.  The  corresponding  projection  of  the 
eigenvector  predicted  by  the  model  onto  the  quantum  mechan¬ 
ical  eigenvector  is  0.94,  indicating  also  a  very  good  prediction 
of  the  vibrational  motion  associated  to  this  mode.  However, 
when  this  potential  was  tested  in  NPT-MD  simulations  in  the 
regime  of  low  temperatures  and  pressures  it  was  found  that  while 
the  unit  cell  size  and  shape  of  nitromethane  were  reasonably 
predicted,  the  methyl  groups  of  the  molecules  in  the  unit  cell 
were  rotated  by  approximately  30°  relative  to  the  experimental 
orientation.  This  fact  indicated  the  need  to  increase  the  barrier 
height  of  the  torsional  H-C-N-O  potentials  in  order  to 
reproduce  accurately  the  crystallographic  structure. 

The  second  set  of  potentials  we  have  considered  is  similar 
to  the  Mgas  model  but  we  have  changed  the  HCNO  torsional 
potential  to  a  3-fold  form  and  increased  the  barrier  parameter 
V<d(H— C-N-O)  to  0.27  kJ/mol,  which  corresponds  to  an 
overall  barrier  for  methyl  rotation  similar  to  that  given  in  ref 
14.  The  complete  list  of  intramolecular  parameters  for  this 
potential  set  is  that  indicated  in  Table  1.  The  modifications  of 
the  corresponding  vibrational  frequencies  for  this  model  are 
shown  in  columns  denoted  as  Ml  and  PI  in  Table  2.  It  can  be 
seen  that  the  most  important  variations  of  the  vibrational 
frequencies  take  place  for  the  mode  vTJ  which  corresponds  to 
the  internal  rotation  of  the  methyl  group.  Overall,  excepting 
modes  vr  and  v3  (which  describes  the  inversion  of  the  methyl 
group),  there  is  good  agreement  between  the  predicted,  the  ab 
initio  calculated,  and  the  experimental24  sets  of  vibrational 
frequencies  with  rms  deviations  of  6.9  and  10.8  cm-1,  respec¬ 
tively.  In  addition,  the  results  indicate  that  the  model  reproduces 
the  quantum  mechanical  eigenvectors  very  well,  with  the 
smallest  projection  values  being  0.92  for  modes  and  vs.  As 
will  be  described  more  fully  hereafter,  MP  calculations  and 
NPT-MD  simulations  using  this  potential  model  reasonably 
reproduce  the  experimental  crystal  structure  and  molecular 
conformations  for  the  entire  range  of  temperatures  and  pressures 
investigated.  Consequently,  all  the  results  presented  in  the  next 
sections  correspond  to  this  choice  of  potential  parameters 
denoted  as  model  Ml. 

A  final  choice  of  torsional  potential  parameters  we  have  tested 
corresponds  to  an  even  larger  value  for  V<j,(H— C— N— O)  of 
0.49  kJ/mol.  Our  results  indicate  that  this  potential  change  has 
only  a  very  small  influence  on  the  predicted  crystallographic 
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TABLE  3:  Lattice  Parameters  and  Energies  Obtained  in  Crystal  Packing  without  Symmetry  Constraints  (the  percentage 
change  in  lattice  and  molecular  parameters  after  energy  minimization  is  determined  as  a  function  of  the  experimental  geometry 
and  is  given  in  parentheses) 


lattice  energy*2 

lattice  parameters^ 

Qc  NB  ES 

a 

b 

c 

a 

p 

Y 

EXP* 

-52.3* 

5.1832 

6.2357 

8.5181 

90.000 

90.000 

90.000 

LMIN 

20.5  -25.0401  -35.0560 

5.2369 

6.2653 

8.6214 

90.001 

89.989 

89.997 

(1.04) 

(0.48) 

(1.21) 

(0.00) 

(-0.01) 

(0.00) 

a / 

Ay 

Az 

A<!>  (deg)* 

A©  (deg) 

A1^  (deg) 

-0.037 

-0.105 

-0.014 

3.44 

3.02 

6.41 

d  Nonbonded  (NB)  and  electrostatic  (ES)  lattice  energies  per  molecule  in  kJ/mol.  b  Lattice  dimensions  a,  b ,  and  c  in  angstroms,  and  angles  a, 

and  y  in  degrees.  c  Cutoff  parameter  as  described  in  text.  d  Experimental  values  from  ref  7.  *  Estimated  using  the  theoretical  heat  of  sublimation  of 

—47.3  kJ/mol  as  provided  by  Politzer.34  f  Change  in  fractional  coordinates  of  molecular  centroids  for  nitromethane  molecule.  8 

Change  in  Euler 

angles  of  molecular  centroids  for  nitromethane. 


parameters  relative  to  the  Ml  model  and  that  this  influence  is 
particularly  limited  to  the  region  of  low  temperatures.  Conse¬ 
quently  we  will  make  only  limited  references  to  this  potential 
choice  in  the  present  work.  For  completeness,  we  indicate  in 
Table  2  columns  M2  and  P2,  the  predicted  vibrational  frequen¬ 
cies  and  the  corresponding  projections  for  this  potential  model. 
As  can  be  seen  only  the  value  of  the  torsional  frequency  vx  is 
increased  by  about  40  cm”1  relative  to  the  Ml  value  while  the 
projection  of  its  eigenvector  on  the  corresponding  ab  initio 
vector  remains  0.94.  The  large  value  of  this  projection,  close 
to  unity,  indicates  that  even  for  this  high  barrier  value  the  model 
is  representing  correctly  the  corresponding  vibrational  atomic 
motions. 

in.  Computational  Details 

Molecular  Packing  Calculations.  A  first  set  of  calculations 
used  to  test  the  semiempirical  intermolecular  potential  energy 
functions  proposed  here  is  based  on  the  use  of  molecular  packing 
calculations,25’26  in  which  the  lattice  energy  of  the  crystal  is 
minimized  with  respect  to  the  structural  degrees  of  freedom  of 
the  crystal.  These  calculations  have  been  done  using  the 
algorithm  proposed  by  Gibson  and  Scheraga27  for  efficient 
minimization  of  the  energy  of  a  fully  variable  lattice  composed 
of  rigid  molecules  and  implemented  in  the  program  LM3N.28 
The  nonbonded  interactions  were  cut  off  with  a  cubic  spline 
function  from  Po  to  Qo ,  to  ensure  the  continuity  of  the  function 
and  its  first  derivative.  Here  o  is  the  value  of  r  in  eq  2  at  which 
V<xj}(r)  =  0  and  dV^rydr  <  0.  The  parameters  P  and  Q ,  which 
specify  the  start  and  the  end  of  the  cubic  feather  (see  refs  1  and 
27  for  details),  were  set  to  20.5  and  20.0,  respectively.  The 
coulombic  potential  terms  of  the  form  given  in  eq  3  are  summed 
over  the  lattice  using  the  Ewald  technique  as  previously 
described.1  Finally,  the  effect  of  pressure  on  the  crystallographic 
parameters  has  been  simulated  by  adding  a  potential  term  of 
the  form  P(V  -  Vo),27  where  Vo  is  the  volume  of  a  suitably 
chosen  unit  cell  at  zero  pressure. 

Constant-Pressure  and  -Temperature  Molecular  Dynam¬ 
ics  Calculations.  The  dynamics  of  nitromethane  as  a  function 
of  temperature  and  pressure  have  been  investigated  using  NPT- 
MD  simulations  based  on  the  total  potential  described  in  eqs 
1—7.  For  comparison  we  have  also  performed  calculations  using 
the  rigid-body  approximation1”6  of  the  molecules  in  the  system. 
In  all  simulations  we  have  used  the  Nose-Hoover  thermo¬ 
stat—  barostat  algorithm29  as  implemented  in  the  program 
DL_POLY_2.0,30  to  simulate  the  crystals  at  various  temper¬ 
atures  and  pressures.  In  this  case,  the  equations  of  motion  for 
molecules  and  the  simulation  cell  are  integrated  using  the  Verlet 


leapfrog  scheme.31  In  the  case  of  rigid-molecules  simulations 
the  molecular  rotational  motion  is  handled  using  Fincham’s 
implicit  quaternion  algorithm.32 

The  MD  simulation  cells  consist  of  boxes  containing  5x4 
x  3  crystallographic  unit  cells.  This  choice  of  the  simulation 
box  ensures  the  use  of  a  cutoff  distance  for  the  intermolecular 
potentials  of  about  10  A.  The  initial  configuration  corresponding 
to  the  lowest  temperature  was  chosen  to  be  identical  to  that  for 
the  low-temperature  experimental  structure.  The  system  was  then 
equilibrated  at  that  temperature  and  atmospheric  pressure.  In 
all  production  runs  done  using  the  Nos6-Hoover  implementation 
for  the  NPT  ensemble,  the  system  was  integrated  for  34000 
time  steps  (1  time  step  =  .75  x  10“15  s),  of  which  4000  steps 
were  equilibration.  In  the  equilibration  period,  the  velocities 
were  scaled  after  every  5  steps  so  that  the  internal  temperature 
of  the  crystal  mimicked  the  imposed  external  temperature.  Then, 
properties  were  calculated  and  accumulated  for  averaging  over 
the  next  30000  integration  steps  in  the  simulation.  In  subsequent 
runs,  performed  at  successively  higher  temperatures  or  pressures, 
the  initial  configurations  of  the  molecular  positions  and  veloci¬ 
ties  were  taken  from  the  previous  simulation  at  the  end  of  the 
production  run. 

The  lattice  sums  were  calculated  subject  to  the  use  of 
minimum-image  periodic  boundary  conditions  in  all  dimen¬ 
sions.31  The  interactions  were  determined  between  the  sites 
(atoms)  in  the  simulation  box  and  the  nearest-image  sites  within 
the  cutoff  distance.  In  these  calculations,  the  coulombic  long- 
range  interactions  were  handled  using  Ewald’s  method.31 

The  main  quantities  obtained  from  these  simulations  were 
the  average  lattice  dimensions  and  the  corresponding  volume 
of  the  unit  cell.  Additional  information  about  the  structure  of 
the  crystal  has  been  obtained  by  calculating  the  radial  distribu¬ 
tion  functions  (RDF)  between  different  atomic  sites.  Such 
quantities  have  been  calculated  from  recordings  done  at  every 
10th  step  during  the  trajectory  integrations. 

TV.  Results  and  Discussions 

A.  Molecular  Packing  Calculations.  The  results  of  MP 
calculations  without  symmetry  constraints  are  presented  in  Table 
3.  As  can  be  seen  the  predicted  structural  lattice  parameters  for 
nitromethane  differ  by  less  than  1.21%  from  the  experimental 
values  reported  by  Trevino  et  al.7  Also,  there  are  small 
translations  but  slightly  larger  rotations  (with  a  maximum 
deviation  of  6.4°)  of  the  molecule  in  the  asymmetric  unit  cell. 
Using  the  HF  set  of  charges,  the  total  lattice  energy  is  predicted 
to  be  -60. 1  kJ/mol.  This  value  compares  acceptably  well  with 
the  lattice  energy  of  -52.3  kJ/mol  that  has  been  estimated  using 
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TABLE  4:  Lattice  Parameters  Obtained  in  NPT-MD  Calculations  for  Nitromethane  as  a  Function  of  Temperature;  the 
Calculated  Thermal  Expansion  Coefficients  (%)  at  250  K  Are  Also  Indicated" 


lattice  dimensions 


r(K) 

b(k) 

c(A) 

a(deg) 

/S(deg) 

y(deg) 

volume(A3) 

4.2* 

5.1832 

6.2357 

8.5181 

275.3125 

78.0* 

5.1983 

6.2457 

8.5640 

278.0476 

228.0* 

5.2440 

6.3200 

8.7300 

289.3303 

4.2 

5.2116(0.54) 

6.3416  (1.69) 

8.6464  (1.50) 

89.998 

89.999 

90.000 

285.7666  (3.79) 

25.0 

5.2134 

6.3608 

8.6935 

90.012 

90.018 

89.996 

288.3000 

50.0 

5.2184 

6.3794 

8.7266 

90.001 

90.034 

90.000 

290.5000 

78.0 

5.2195  (0.40) 

6.4032  (2.52) 

8.7543  (2.22) 

89.991 

89.996 

89.991 

292.5833  (5.22) 

100.0 

5.2222 

6.4173 

8.7734 

90.001 

89.962 

90.007 

294.0166 

125.0 

5.2267 

6.4329 

8.7979 

89.987 

89.994 

90.002 

295.8166 

150.0 

5.2298 

6.4605 

8.8204 

89.974 

89.991 

90.010 

298.0166 

175.0 

5.2373 

6.4889 

8.8606 

89.995 

90.005 

89.987 

301.1166 

200.0 

5.2452 

6.5167 

8.8808 

90.016 

90.012 

90.020 

305.3166 

228.0 

5.2534  (0.18) 

6.5481  (3.60) 

8.9195  (2.17) 

90.015 

90.020  ’ 

90.043 

303.5667  (4.92) 

250.0 

5.2620 

6.5691 

8.9561 

90.001 

89.993 

89.946 

306.8167 

X 

68.9  x  10“6 

181.7  x  10"6 

131.9  x  10"6 

380.5  x  10-6 

a  The  values  in  parentheses  represent  the  percentage  differences  relative  to  the  available  experimental  results.  *  Experimental  data  from  ref  7. 
c  The  units  for  the  linear  and  volume  expansion  coefficients  are  K'1. 


the  relation  f/lat  «  -Atfsubi  -  2 RTP  In  this  analysis  we  have 
used  as  the  heat  of  sublimation  for  nitromethane  the  value  of 
47.2  kJ/mol  estimated  theoretically  by  Politzer.34 

B.  NPT-MD  Calculations.  NPT-MD  simulations  of  nitro¬ 
methane  were  used  to  predict  the  crystal  structure  of  ni¬ 
tromethane  over  a  large  range  of  temperatures  and  pressures. 
We  will  first  examine  thermal  effects  on  the  crystallographic 
cell  and  molecular  parameters  at  1  atm,  and  compare  them 
against  experimental  values. 

Bl.  Temperature  Effects .  Average  unit  cell  edge  lengths  and 
volumes  determined  from  NPT-MD  simulations  for  the  flexible 
model  for  temperatures  ranging  from  4.2  to  250  K  are  given  in 
Table  4.  These  data  are  represented  in  Figure  2,  along  with  the 
corresponding  experimental  values  and  the  NPT-MD  results 
obtained  using  the  rigid-body  approximation.  For  the  flexible 
model,  the  NPT-MD  lattice  dimensions  obtained  at  T  =  4.2  K 
are  in  very  close  agreement  with  those  determined  in  the  MP 
calculations  and  with  the  experimental  values.7  At  this  temper¬ 
ature,  the  percentage  differences  between  the  predicted  and  the 
experimental  lattice  dimensions  are  0.54%,  1.69%,  and  1.50% 
for  the  a,  b,  and  c  axes,  respectively.  The  deviation  of  the 
predicted  unit  cell  volume  from  the  experimental  value  at  this 
temperature  is  3.79%.  This  good  agreement  is  decreased  only 
slightly  with  an  increase  in  temperature.  At  T  =  78  K,  the  cell 
edge  lengths  and  the  unit  cell  volume  agree  with  the  experi¬ 
mental  values  to  within  2.52%  and  5.22%,  respectively.  At  the 
highest  temperature  where  experimental  data  are  available  (228 
K),  the  deviations  from  experiment  for  the  cell  edges  are  less 
than  3.60%  while  the  difference  is  4.92%  for  the  unit  cell 
volume.  Also,  the  unit  cell  angles  remain  close  to  90.0°  for  the 
entire  temperature  range  analyzed,  as  expected  for  the  orthor¬ 
hombic  symmetry  of  the  crystal.  For  comparison,  calculations 
using  the  M2  potential  model  indicate  slightly  larger  deviations 
for  the  predicted  lattice  parameters  from  experimental  values 
of  1.31%,  1.03%,  and  1.72%  at  7=4.2K  and  0.60%,  3.61%, 
and  2.13%  at  T  =  228  K  for  the  a,  b,  and  c  axes,  respectively. 

The  linear  and  volume  thermal  expansion  coefficients  ex¬ 
tracted  from  the  data  are  also  given  in  Table  4.  The  expansion 
of  the  lattice  is  highly  anisotropic  with  the  largest  length  changes 
along  the  b  and  c  axes.  These  effects  can  be  partially  attributed 
to  the  internal  rotation  of  the  methyl  group.  To  illustrate  our 
argument,  consider  the  limiting  case  in  which  a  nitromethane 
molecule  is  arranged  in  a  unit  cell  such  that  the  C— N  bond  is 
parallel  to  the  a  axis.  Changes  in  the  positions  of  the  hydrogen 
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Figure  2.  Variation  of  the  lattice  dimensions  (a)  and  unit  cell  volume 
(b)  with  temperature  from  NPT-MD  simulations  using  the  rigid-body 
approximation  (MD,  R)  and  the  flexible  model  (MD,  F).  The  available 
experimental  data  from  ref  7  (EXP)  are  also  indicated. 

atoms  due  to  the  rotational  motion  of  the  methyl  group  about 
the  C— N  bond  would  occur  in  the  b— c  plane  only.  Thus,  the 
increasingly  large  librational  motion  occurring  with  increases 
in  temperature  would  affect  only  the  cell  dimensions  associated 
with  the  b  and  c  axes.  In  the  actual  crystal,  the  arrangement  of 
the  molecules  in  the  unit  cell  is  such  that  the  C— N  bonds  of  all 
four  are  nearly  perpendicular  to  the  c  axis  (as  can  be  seen  in 
Figure  1).  The  angles  that  each  C-N  bond  makes  with  the  axes 
of  the  crystal  cell  are  40,  51,  and  86°  for  the  a,  b,  and  c  axes, 
respectively.  This  alignment  of  the  molecules  relative  to  these 
axes  indicates  that  rotations  of  the  methyl  group  about  the  C— N 
bond  would  more  strongly  affect  the  b  and  c  axes.  Indeed,  by 
comparing  the  data  in  Figure  2a,  it  is  clear  that  the  major 
differences  between  NPT-MD  data  from  the  rigid  and  flexible 
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TABLE  5:  Comparison  of  the  Average  Orientational  Parameters  (fractional  coordinates  sx,  sy,  sz  and  Euler  angles  0,  «I>,  V) 
for  the  Four  Molecules  in  the  Unit  Cell  and  of  the  Internal  Geometrical  Parameters  of  Nitromethane  Molecule  as  Obtained 
from  Trajectory  Calculations  at  Different  Temperatures;  Corresponding  Experimental  Values  (where  available)  Are  Also 
Indicated0 


parameter 

temperature(K) 

4.2 

expf7^ 

4.2 

NPT-MD 

4.2 

ideal 

78 

expt^ 

78 

NPT-MD 

78 

ideal 

150 

NPT-MD 

150 

ideal 

sxl 

0.3493 

0.3513 

0.3475 

0.3483 

0.3484 

syl 

0.4196 

0.4284 

0.4188 

0.4288 

0.4290 

szl 

0.3750 

0.3682 

0.3803 

0.3722 

0.3740 

sx2 

0.1507 

0.1487 

0.1487 

0.1525 

0.1509 

0.1517 

0.1512 

0.1515 

sy2 

-0.4196 

-0.4284 

-0.4284 

-0.4188 

-0.4291 

-0.4288 

-0.4290 

-0.4290 

sz2 

-0.1250 

-0.1318 

-0.1318 

-0.1197 

-0.1279 

-0.1278 

-0.1262 

-0.1260 

sx3 

-0.3493 

-0.3513 

-0.3513 

-0.3475 

-0.3494 

-0.3483 

-0.3488 

-0.3485 

sy3 

-0.0804 

-0.0717 

-0.0716 

-0.0812 

-0.0711 

-0.0712 

-0.0705 

-0.0710 

sz3 

0.1250 

0.1317 

0.1318 

0.1197 

0.1277 

0.1278 

0.1260 

0.1260 

sx4 

-0.1507 

-0.1488 

-0.1487 

-0.1525 

-0.1513 

-0.1517 

-0.1516 

-0.1515 

sy4 

0.0804 

0.0716 

0.0716 

0.0812 

0.0708 

0.0712 

0.0703 

0.0710 

sz4 

-0.3750 

-0.3683 

-0.3682 

-0.3803 

-0.3719 

-0.3722 

-0.3739 

-0.3740 

©i 

90.8 

94.4 

96.6 

92.9 

92.2 

4>1 

53.5 

51.9 

51.6 

53.2 

53.7 

MU 

238.9 

245.3 

240.1 

245.4 

245.4 

02 

-90.8 

-94.4 

-94.4 

-96.6 

-92.9 

-92.9 

-92.4 

-92.2 

4>2 

53.5 

52.0 

51.9 

51.6 

53.1 

53.2 

53.7 

53.7 

W2 

238.9 

245.3 

245.3 

240.1 

245.5 

245.4 

245.4 

245.4 

©3 

89.2 

85.5 

85.6 

83.4 

87.2 

87.1 

87.7 

87.8 

4>3 

-53.5 

-52.0 

-51.9 

-51.6 

-53.1 

-53.2 

-53.8 

-53.7 

58.9 

65.3 

65.3 

60.1 

65.4 

65.4 

65.4 

65.4 

©4 

-89.2 

-85.6 

-85.6 

-83.4 

-87.1 

-87.1 

-87.7 

-87.8 

4>4 

-53.5 

-52.0 

-51.9 

-51.6 

-53.2 

-53.2 

-53.7 

-53.7 

\p4 

58.9 

65.3 

65.3 

60.1 

65.4 

65.4 

65.4 

65.4 

temperature(K) 

4.2 

4.2 

78 

78 

150 

parameter 

expt* 

NPT-MD 

expt.* 

NPT-MD 

NPT-MD 

C-N 

1.481 

1.497 

1.488 

1.499 

1.499 

C— H 

1.098 

1.088 

1.077 

1.089 

1.090 

n-o3 

1.209 

1.230 

1.198 

1.230 

1.231 

n-o4 

1.223 

1.229 

1.226 

1.230 

1.231 

HCH 

111.6 

110.7 

111.2 

110.6 

110.5 

HCN 

107.3 

108.2 

107.7 

108.2 

108.2 

cno3 

118.9 

118.1 

119.1 

117.8 

117.8 

cno4 

117.8 

117.7 

116.8 

117.8 

117.9 

O3NO4 

123.3 

124.3 

124.1 

124.2 

124.2 

0  In  the  case  of  molecules  2-4  the  ideal  orientational  parameters  obtained  based  on  P2i2]2t  space  group  symmetry  operators  and  the  orientation 
parameters  of  molecule  1,  positioned  in  the  asymmetric  unit,  are  given.  b  Experimental  values  from  ref  7.  c  Unit  cell  fractional  values  have  been 
shifted  to  fall  within  the  range  from  - 1/2  to  1/2. 


potential  models  are  in  the  expansions  along  the  b  and  c  axes 
while  the  variations  for  the  case  of  the  a  axis  are  practically 
absent.  The  additional  expansion  due  to  the  methyl  group 
rotation  can  be  also  seen  in  Figure  2b  where  the  volume  values 
obtained  using  the  flexible  model  are  consistently  larger  than 
the  corresponding  data  obtained  for  the  rigid  model. 

In  Table  5  we  compare  first  the  results  of  the  average 
fractional  coordinates  and  orientational  Euler  parameters  of  the 
four  molecules  in  the  unit  cell  with  the  corresponding  experi¬ 
mental  data  (where  available).  These  values  are  averaged  over 
time  and  all  unit  cells  in  the  simulation  box.  Additionally,  in 
the  same  table  we  provide  the  “ideal”  orientational  Euler 
parameters  for  three  of  the  four  molecules  of  the  unit  cell  relative 
to  the  orientation  of  the  reference  molecule  in  the  asymmetric 
unit  cell.  These  parameters  have  been  calculated  assuming  the 
symmetry  operations  of  the  P2i2i2i  space  group.  In  this  way 
we  can  analyze  the  degree  of  deviation  of  the  predicted  crystal 
structure  ffom  the  P2i2i2i  symmetry. 

It  is  evident  that  increasing  the  temperature  ffom  4.2  to  150 
K  does  not  produce  any  significant  displacement  of  the 
molecular  centers  of  mass  and  that  the  degree  of  rotational 
disorder  is  small.  The  largest  deviations  of  the  system  from  the 


orientational  parameters  corresponding  to  a  perfect  crystal  with 
P2i2j2i  symmetry  are  0.0011  and  0.2°  for  the  fractional 
coordinates  and  Euler  angles,  respectively.  A  similar  conclusion 
is  obtained  for  temperatures  between  150  and  250  K  (not 
shown).  At  T  —  4.2  K,  the  largest  deviation  between  the 
experimental  and  predicted  molecular  orientations  is  about  6.4° 
for  the  Euler  angle  W  (see  Table  5),  in  agreement  with  previous 
MP  findings.  At  T  =  78  K  this  maximum  deviation  decreases 
to  about  5.4°. 

Additional  support  for  the  small  degree  of  translation  of  the 
molecules  inside  the  unit  cell  with  the  temperature  increase  can 
be  obtained  ffom  the  C—C  RDFs  given  in  Figure  3a.  As  can 
be  seen,  the  RDFs  at  these  temperatures  correspond  to  well- 
ordered  structures  with  correlation  at  long  distances.  The 
positions  of  the  major  peaks  do  not  change  significantly  and 
the  main  temperature  effect  is  the  broadening  of  the  peaks  with 
partial  overlapping  of  some  of  them. 

The  good  agreement  of  the  structural  molecular  parameters 
predicted  by  the  model  with  experiment  can  be  also  observed 
by  analyzing  the  internal  coordinates  (bond  lengths  and  angles) 
indicated  in  Table  5.  At  4.2  K  the  maximum  deviations  for  bond 
lengths  is  about  1.7%  and  about  1.0%  for  bond  angles.  This 
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Figure  3.  Radial  distribution  functions  for  C***C  and  0-*H  as  functions 
(c-d)atr=  293  K. 

level  of  agreement  is  particularly  good  taking  into  account  that 
in  our  potential  parametrization,  we  have  considered  the  internal 
geometrical  parameters  of  the  nitromethane  molecule  determined 
from  gas-phase  ab  initio  calculated  values  rather  than  parameters 
corresponding  to  nitromethane  in  the  crystalline  phase.  These 
results  suggest  also  a  small  influence  of  the  lattice  field  on  the 
internal  parameters  of  the  nitromethane  molecule.  Changes  in 
predicted  molecular  structural  parameters  with  temperature  are 
negligible,  again  in  agreement  with  experiment.  The  most 
significant  change  appears  due  to  the  internal  rotation  of  the 
methyl  group  as  a  function  of  temperature.  This  is  evident  by 
inspecting  the  behavior  of  the  radial  distribution  functions 
(RDFs)  for  the  O—H  intermolecular  bonds  given  in  Figure  3b. 
There  is  a  rapid  destruction  of  the  correlation  at  long  distances 
with  the  increase  of  temperature  and  consequently  the  indi¬ 
viduality  of  the  H  atoms  in  the  C— H—O  bonds  is  lost.  These 
effects  can  be  understood  as  being  due  to  the  increase  in 
rotational  disorder  of  the  methyl  group  with  the  increase  of 
temperature. 

TTie  increase  in  the  librational  motion  of  the  methyl  group 
with  temperature  is  illustrated  in  Figure  4.  This  figure  shows 
the  cumulative  distributions  of  the  H;— Q—  N2— O3  (/  =  5—7) 
dihedral  angles  at  temperatures  ranging  from  4.2  to  150  K.  These 
distributions  have  been  determined  from  the  molecular  con¬ 
figurations  recorded  at  every  tenth  step  during  trajectories  of 
16000  steps  (12  ps).  In  this  figure,  the  distributions  are  peaked 
at  the  angles  corresponding  to  the  equilibrium  orientation  of 
the  methyl  group  in  the  low-temperature  structure,  namely, 


r(A) 


r(A) 

of  temperature  (a— b)  at  atmospheric  pressure  and  as  function  of  pressure 

-78°,  42°,  and  161°  respectively.  These  angles  are  within  10° 
from  those  corresponding  to  the  experimental  structure  deter¬ 
mined  by  the  neutron  diffraction  technique.7  As  can  be  seen  in 
Figure  4  the  increase  of  temperature  does  not  change  the  position 
of  these  peaks.  However,  the  distributions  broaden  with  tem¬ 
perature,  indicating  an  increase  in  librational  motion  with 
temperature.  For  the  case  of  the  potential  model  with  the 
increased  HCNO  torsional  barrier  (M2)  we  have  found  that  the 
agreement  between  the  predicted  (—85°,  34°,  154°)  and 
experimental  (-88°,  31°,  151°)  torsional  angles  is  even  better, 
with  deviations  within  3°.  Similarly,  in  this  case  the  increase 
of  temperature  did  not  modify  the  position  of  the  peaks  in  the 
cumulative  distributions. 

The  magnitudes  of  the  methyl  group  libration  have  been 
previously  determined  experimentally  for  deuterated  ni¬ 
tromethane  in  the  temperature  range  4.2—125  K  using  neutron 
diffraction  powder  experiments.7  As  can  be  seen  in  Figure  5 
where  the  experimental  values  are  compared  with  our  results 
the  agreement  is  quite  good  for  the  entire  temperature  range 
25-125  K.  For  example,  the  experimental  results  indicate  that 
at  25  and  125  K  the  rms  amplitudes  of  libration  of  the  CD3 
group  are  ~13°  and  25°,  respectively.  Our  NPT-MD  calculations 
of  deuterated  nitromethane  predict  rms  amplitudes  of  libration 
of  12°  at  25  K  and  23°  at  125  K  (see  Figure  5).  Comparison  of 
the  amplitude  of  libration  at  4.2  K  with  experiment  is  not 
appropriate,  since  quantum  mechanical  effects  would  be  sig¬ 
nificant  at  this  temperature.  The  experiments7  suggest  that  the 
amplitude  of  libration  at  this  temperature  is  as  large  as  that  at 
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Figure  4.  Distribution  of  the  H, — C\—  N2— O3  (1  —  5,6,7)  dihedral 
angles  for  all  nitromethane  molecules  in  the  simulation  box  as  a  function 
of  temperature  (a)  and  pressure  (b). 


Figure  5.  Root-mean-square  amplitude  of  libration  of  the  methyl  group 
in  CD3NO2  as  a  function  of  temperature. 

25  K;  our  classical  simulations  predict  a  value  that  is  smaller 
by  ~7°.  Previous  quasielastic  and  inelastic  neutron  scattering 
studies8  performed  on  protonated  nitromethane  have  provided 
additional  information  about  the  temperature  dependence  of  the 
rotational  reorientation  of  the  methyl  groups.  By  assuming  a 
3-folded  reorientation  model  and  an  Arrhenius  dependence 
(^eE* /kT)  0f  the  mean  residence  times  of  methyl  groups  on  the 
temperature  7,  an  activation  barrier  of  234  cal/mol  to  the  methyl 
group  rotation  has  been  determined.  However,  more  recent 
experimental  studies  based  on  infrared  and  Raman  measure¬ 
ments  have  indicated  higher  activation  energies  with  values  of 
401  cal/mol16  and  576  cal/mol,10  respectively.  We  monitored 


i/T  (ic1) 


Figure  6.  The  mean  residence  times  r  of  the  methyl  group  versus 
inverse  temperature.  The  activation  energy  of  387  cal/mol  is  obtained 
from  the  slope  of  the  least-squares  linear  fit  to  the  data  in  the  region  of 
higher  temperatures  (125-250  K)  indicated  by  a  solid  line. 

the  frequency  of  120°  methyl  group  rotations  for  the  duration 
of  12  ps  in  trajectories  integrated  at  temperatures  ranging  from 
4.2  to  250  K.  From  these  trajectories  we  have  observed  that 
methyl  groups  undergo  large  vibrational  amplitudes,  indicative 
of  small  barrier  heights  for  torsional  motions  in  agreement  with 
experimental  findings. 

To  obtain  a  more  quantitative  description  of  the  internal 
rotational  motion  we  have  considered  the  evaluation  of  tem¬ 
perature  effects  on  the  rates  of  reorientation  of  the  methyl  group. 
A  “rotational  jump”  was  defined  in  the  following  manner.  At 
the  beginning  of  each  trajectory,  the  H/— Ci—  N2—O3  angles  (i 
=  5-7)  of  each  molecule  were  calculated.  Each  angle  was 
assigned  to  one  of  the  three  minima  in  the  3-fold  rotational 
potential  by  virtue  of  its  proximity  to  the  minima,  and  the  time 
denoted  as  to.  The  minima  corresponded  to  the  values  of  the 
H/-C1-N2-O3  angles  (1  =  5—7)  in  the  low-temperature 
structure,  i.e.,  —78°,  42°,  and  161°,  respectively.  The  assignment 
of  these  angles  to  the  minima  are  made  at  each  subsequent  step 
in  the  trajectory  and  compared  to  the  original  assignment  at  r0- 
When  the  assignments  for  all  three  of  the  H— Q— N2-O3  angles 
differ  from  the  original  assignments  at  to,  it  is  assumed  that  a 
rotational  jump  has  occurred.  At  this  point,  the  time  interval 
between  the  time  of  the  new  assignments  and  to  was  recorded 
The  value  of  to  was  reset  to  the  time  of  the  new  assignments, 
and  the  process  repeated  for  the  duration  of  the  trajectory. 
Methyl  group  reorientation  rates  at  each  temperature  were 
determined  by  extracting  first-order  rate  coefficients  from  linear 
fits  of  plots  of  ln(P)  versus  time,  where  P  is  the  fraction  of 
nitromethane  molecules  that  have  not  jumped  at  time  /.  The 
lifetimes  indicate  first-order  behavior.  An  Arrhenius  plot  of  the 
residence  times,  which  are  the  inverse  of  the  reorientation  rates, 
is  shown  in  Figure  6.  Fitting  a  line  to  all  the  points  in  Figure  6 
predicts  Ez  =  241  cal/mol.  However,  the  nonlinearity  of  the 
Arrhenius  plot  is  a  result  of  the  statistical  errors  in  the  low- 
temperature  rates,  and  if  we  base  our  prediction  on  the  7  points 
for  the  highest  temperaures  (100—250  K)  we  predict  Ea  =  387 
cal/mol.  These  values  are  in  accord  with  the  reported  experi¬ 
mental  activation  energies  (234—576  cal/mol).8, iai6 

B2.  Pressure  Effects.  The  calculated  lattice  dimensions  at 
different  pressures  are  given  in  Table  6  and  a  visual  comparison 
with  the  experimental  data  is  presented  in  Figure  7.  For  the 
entire  pressure  range  investigated  (0.3—7  GPa)  the  differences 
between  the  calculated  NPT-MD  data  using  the  flexible  potential 
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TABLE  6:  Lattice  Parameters  Obtained  in  Crystal  Packing  and  NPT-MD  Calculations  for  Nitromethane  as  a  Function  of 


Pressure 


press  (GPa) 

NPT-MD 

molecular  packing 

lattice  lengths  (A) 

unit  cell 
volume  (A3) 

lattice  lengths  (A) 

unit  cell 
volume  (A3) 

a 

b 

c 

a 

b 

c 

0.30 

5.2070 

6.4829 

8.8597 

299.0666 

5.2068 

6.1989 

8.5628 

276.3758 

0.33 

5.2006 

6.4874 

8.8416 

298.2833 

5.2040 

6.1929 

8.5576 

275.7927 

0.50 

5.1713 

6.4236 

8.7800 

291.6500 

5.1861 

6.1561 

8.5172 

271.7649 

1.00 

5.0997 

6.3001 

8.6595 

278.2166 

5.1397 

6.0786 

8.4295 

262.9498 

2.00 

4.9827 

6.1706 

8.4936 

261.1500 

5.0656 

5.9716 

8.3055 

250.5347 

3.00 

4.8963 

6.0797 

8.3826 

249.5166 

5.0051 

5.8988 

8.2150 

241.6301 

3.50 

4.8603 

6.0448 

8.3392 

245.0000 

4.9792 

5.8684 

8.1772 

237.9397 

4.00 

4.8245 

6.0139 

8.2955 

240.6833 

4.9536 

5.8447 

8.1421 

234.6472 

5.00 

4.7652 

5.9600 

8.2265 

233.6500 

4.9054 

5.8085 

8.0780 

228.8875 

5.45 

4.7426 

5.9412 

8.2001 

231.0500 

4.8843 

5.7969 

8.0517 

226.5996 

6.00 

4.7116 

5.9145 

8.1666 

227.5833 

4.8527 

5.7982 

8.0161 

223.9448 

6.50 

4.6898 

5.8925 

8.1404 

224.9666 

4.7698 

5;9135 

7.9484 

221.2719 

7.00 

4.6666 

5.8733 

8.1162 

222.4500 

4.7171 

5.9801 

7.8992 

218.8226 

Figure  7.  Pressure  variation  of  the  lattice  dimensions  (a),  unit  cell 
volume  (b),  and  volume  compression  V/Vo  on  the  external  pressure. 
The  calculated  data  from  molecular  packing  (MP),  NPT-MD  simulations 
using  the  rigid-body  approximation  (MD,R)  and  the  flexible  model 
(MD,F)  are  represented.  The  corresponding  experimental  data  from  refs 
9  (EXP)  and  12  (EXP2),  respectively,  are  also  indicated. 


and  the  experimental  crystallographic  parameters  are  small.  By 
considering  as  reference  the  experimental  data  determined  by 
Cromer  et  al.9  we  find  that  in  the  low-pressure  region  (0.3  GPa) 
the  percentage  errors  for  lattice  dimensions  a,  b,  and  c  are  0.3%, 
3.4%,  and  2.9%,  respectively.  A  similar  good  correspondence 
is  found  in  the  high-pressure  region  (6.0  GPa)  where  the 
percentage  errors  are  0.4%,  3.5%,  and  2.4%  for  the  a,  b,  and  c 
lattice  dimensions,  respectively.  Moreover,  for  the  entire  pres¬ 
sure  range  investigated  the  NPT-MD  predicted  volume  com¬ 
pression  closely  follows  the  corresponding  variation  observed 


experimentally  (see  Figure  7b).  In  Figure  7b  we  have  also 
represented  the  experimental  values  reported  by  Yarger  and 
Olinger.12  As  can  be  seen  in  this  case  the  agreement  with  our 
calculated  unit  cell  volumes  is  even  better  than  that  found  with 
data  provided  by  Cromer  et  al.,9  with  decreasing  deviations  from 
about  3.6%  at  1  GPa  to  1.4%  at  7  GPa. 

The  MP  results  obtained  using  the  rigid  body  approximation 
generally  follow  the  predictions  obtained  using  the  flexible 
potentials,  particularly  in  the  low-pressure  region.  Also,  the 
agreement  between  the  NPT-MD  data  obtained  using  the  rigid 
molecular  assumption  is  surprisingly  good  with  NPT-MD  data 
using  the  flexible  model.  This  fact  indicates  that  for  the  pressure 
region  investigated,  the  compression  of  the  lattice  is  almost 
entirely  due  to  a  reduction  of  intermolecular  distances.  These 
effects  can  be  also  seen  from  the  inspection  of  the  normalized 
dependence  of  the  unit  cell  volume  on  pressure  given  in  Figure 
7c.  In  this  plot  the  experimental  curve9  and  those  predicted  by 
the  flexible  and  rigid  models  are  practically  superimposed  for 
the  entire  range  of  pressures  investigated.  The  values  obtained 
in  the  MP  calculations  start  to  deviate  more  significantly  from 
the  experimental  results,  particularly  for  pressures  above  5  GPa. 
These  facts  indicate  that  the  lattice  compressibility  becomes  less 
well  represented  in  these  calculations  as  the  pressure  is 
increased.  However,  the  present  set  of  results  support  our 
previous  findings6  that  in  the  region  of  low  to  moderate  pressures 
(~5  GPa)  the  MP  calculations  can  be  used  as  an  alternative 
tool  to  describe  the  changes  of  the  unit  cell  geometrical 
parameters  but  at  a  fraction  of  the  computational  cost  involved 
in  MD  simulations. 

Further  comparison  of  the  predicted  and  experimental  lattice 
parameters  can  be  obtained  by  extrapolating  to  zero  pressure 
and  293  K  the  results  of  different  studies,  which  in  principle 
should  coincide.  As  in  previous  experimental  studies,7’9’12  the 
extrapolation  is  taken  because  at  room  temperature  nitromethane 
is  in  the  liquid  phase.  The  coefficients  of  cubic  fits  in  pressure 
of  the  predicted  lattice  parameters  and  unit  cell  volume  to  the 
data  predicted  by  our  flexible  potential  are  given  in  Table  7.  In 
Table  8  we  compare  the  zero-pressure  extrapolated  data  obtained 
in  this  work  to  those  of  previous  experimental  investigations. 
Particularly,  we  have  considered  both  the  low-temperature  data 
obtained  by  Trevino  et  al7  and  extrapolated  to  room  temperature 
as  well  as  the  high-pressure  values  determined  previously9’12 
and  extrapolated  to  zero  pressure.  As  can  be  seen  in  Table  8, 
our  crystal  lattice  dimensions  are  within  2.2%,  2.9%,  and  3.4% 
from  the  data  given  in  refs  7,  9,  and,  12,  respectively.  This  level 
of  agreement  is  satisfactory  taking  into  account  the  relative  large 
deviations  between  different  sets  of  experimental  values. 
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TABLE  7:  Coefficients  of  the  Cubic  Fit  in  Pressure  (GPa)  of  the  Lattice  Constants  and  Unit  Cell  Volume 


Sorescu  et  al. 


_ ao 

a( A)  5.246  A 

b( A)  6.537  A 

c(A)  8.911  A 

V(A3)  305.17  A 


_ ai _ 

—  1.601  x  10"1  AGPa*"1 
-2.475  x  10"1  A  GPa-*1 
-2.711  x  10_1  A  GPa-1 
-2.885  x  10+1  A3  GPa-1 


_ 02 _ 

1.704  x  10~2  A  GPa"2 
3.836  x  10-2  A  GPa~2 
3.788  x  10"2  A  GPa"2 
4.123  A3  GPa-2 


_ 03 _ 

-0.858  x  10~3  A  GPa"3 
-2.379  x  10~3  A  GPa-3 
-2.210  x  10"3  A  GPa"3 
-2.428  x  10_1  A3  GPa-3 


TABLE  8:  Comparison  of  the  Lattice  Constants  for  Solid 
Nitromethane  Extrapolated  to  293  K  and  0  GPa  and  of  the 
Murnagham  Equation  Coefficients  (eq  8)  from  Present 
Calculations  and  Various  Studies 


parameter 

ref  7 
(low 

temperature) 

ref  9 
(high 
pressure) 

ref  10 
(high 
pressure) 

present 

work 

ao(A) 

5.270 

5.227 

5.197 

5.238 

i>o(A) 

6.375 

6.330 

6.292 

6.514 

Co(A) 

8.832 

8.758 

8.747 

8.889 

Vo(A3) 

296.72 

289.77 

286.02 

302.80 

vwA3) 

292.7 

291.87 

310.35 

Bo(GPa) 

7.0 

9.14 

6.78 

% 

5.7 

6.04 

5.88 

Another  set  of  parameters  we  have  used  to  assess  the  accuracy 

of  our  potential  is  represented  by  the  coefficients  of  the 
Mumaghan  equation:935 


which  has  been  used  to  fit  the  dependence  of  the  unit  cell 
volume  on  pressure.  In  eq  8,  V  is  the  volume  at  pressure  Py  V0 
is  the  volume  at  P  =  0  (and  allowed  to  be  varied  in  the  fitting 
procedure).  Bo  is  the  bulk  modulus  at  zero  pressure,  and  B'q  — 
dBo/dP.  In  Table  8  we  present  the  best-fit  parameters  obtained 
from  our  NPT-MD  calculated  values  together  with  the  corre¬ 
sponding  experimental  data.  Our  results  clearly  support  the 
experimental  findings  obtained  by  Cromer  et  al.9  with  relative 
percentage  deviations  of  3.2%  and  3.1%  for  the  bulk  modulus 
and  its  pressure  derivative  values,  respectively. 

From  the  compression  data  for  our  lattice  parameters  we  have 
also  determined  the  variation  of  cell  edge  moduli,  -X  3P/9X, 
as  a  function  of  pressure  (see  Figure  8).  Our  results  indicate 
that  the  compressibilities  of  the  b  and  c  axes  are  close  to  each 
other,  the  c  axis  being  the  least  compressible  at  high  pressures. 
These  findings  are  similar  to  those  obtained  by  Cromer  et  al.9 
and  represented  for  comparison  in  Figure  8.  However,  the 
differences  in  compressibilities  for  the  b  and  c  axes  with 
pressure  are  more  pronounced  in  the  experimental  case  than 
our  data  show.  Additionally,  we  find  that  the  a  axis  is  the  least 
compressible  at  low  pressures  but  becomes  the  most  compress¬ 
ible  for  pressures  above  2  GPa.  This  trend  is  also  present  in  the 
Cromer  et  al.  data9  and  has  been  also  noted  previously  by  Yarger 
and  Olinger.12  These  authors12  have  suggested  that  the  behavior 
of  the  compressive  response  of  the  unit  cell,  particularly  with 
regard  to  the  a  axis,  could  be  explained  by  the  formation  of 
“chains”  of  molecules,  linked  by  bonds  between  the  nitro  and 
methyl  groups  of  the  neighboring  molecules  as  they  are 
compressed.  Unfortunately,  there  is  no  additional  theoretical  or 
experimental  support  to  identify  such  a  mechanism  and  further 
studies,  particularly  based  on  ab  initio  calculations,  are  necessary 
to  explore  the  formation  of  such  chainlike  structures.  The  only 
available  data  are  related  to  a  density  functional  theoretical  study 
performed  at  the  B3LYP/6-31G**  level  to  describe  the  interac¬ 
tion  between  two  nitromethane  molecules  separated  by  distances 
ranging  from  1.6  to  3.5  A.36  The  molecules  were  initially 
arranged  with  the  methyl  group  on  one  molecule  facing  the  nitro 


Figure  8.  The  variation  of  the  cell  edge  moduli  as  a  function  of 
pressure  from  the  present  calculation  and  from  experimental  data  given 
in  ref  9. 


group  on  the  second  (an  arrangement  that  would  correspond  to 
compression  along  the  a  axis).  As  the  molecules  were  brought 
closer  together,  the  potential  energy  first  showed  a  steep 
increase,  and  the  methyl  group  was  “flattened”  to  the  point  that 
the  hydrogen  and  carbon  atoms  were  located  in  the  same  plane. 
Further  decreases  in  the  intermolecular  distance  led  to  reaction 
to  form  products  HONO,  CH2O,  and  CH3O.  Thus,  these 
calculations  do  not  support  the  hypothesis  that  “chains”  of 
molecules  are  formed  through  new  bonding  of  intermolecular 
methyl—nitro  group  interactions  upon  compression  of  the  a  axis. 
However,  more  detailed  analyses  corresponding  to  the  specific 
interactions  that  take  place  in  crystals  are  necessary.  The  change 
in  the  ordering  of  the  most  compressible  axes  appears  to  occur 
near  the  pressure  at  which  the  equilibrium  orientation  of  the 
methyl  group  has  changed  to  45°  relative  to  the  crystal  at  low 
pressure.9  This  pressure-induced  rotational  reorientation  provides 
another  test  of  our  interaction  potential,  the  results  of  which 
we  discuss  next. 

Changes  in  the  rotational  motion  of  the  methyl  group  as  a 
function  of  pressure  can  be  observed  through  RDFs  for  the 
O—H  intermolecular  bonds.  We  have  already  seen  from 
previous  sections  that  RDFs  for  0**H  intermolecular  bonds 
indicate  an  increase  of  the  degree  of  rotational  disorder  of  the 
methyl  group  as  the  temperature  increases  (at  ambient  pressure). 
The  RDFs  for  O—H  intermolecular  bonds  follow  a  reverse  trend 
when  the  pressure  is  increased  (at  ambient  temperature),  as 
shown  in  Figure  3d.  In  this  case  the  spectrum  becomes  better 
resolved  with  pressure  increases,  indicating  a  significant  de¬ 
crease  of  the  amplitude  of  methyl  rotations  around  the  equilib¬ 
rium  positions. 

An  illustration  of  the  dynamic  behavior  of  the  methyl  group 
librations  with  increasing  pressure  is  presented  in  Figure  9.  In 
this  figure,  the  time  histories  of  the  H;-Ci~ N2— O3  (i  =  5 —7) 
torsional  angles  of  one  of  the  molecules  in  the  crystal  at  two 
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Figure  9.  The  time  history  of  the  variation  of  H, — Cj— N2— O3  (i  = 
5,6,7)  dihedral  angles  for  one  of  the  nitromethane  molecules  from 
trajectories  performed  at  P  =  0.3  GPa  (a),  and  P  —  7.0  GPa  (b),  and 
T  =  293  K.  For  visualization  purposes  some  of  the  angles  have  been 
shifted  by  ±360°. 

different  pressures  (T  =  293  K)  are  represented-  In  the  low- 
pressure  case,  the  methyl  group  rotation  is  slightly  hindered,  as 
methyl  group  “jumps”  from  one  equilibrium  configuration  to 
another  are  seen  at  times  0.5  and  5  ps  in  the  time  histories.  In 
the  high-pressure  regime  the  torsional  angles  oscillate  around 
well-defined  equilibrium  values  and  the  jumps  between  the 
equilibrium  wells  are  absent.  Also,  the  equilibrium  configuration 
for  the  methyl  group  changes  with  pressure,  as  illustrated  by 
the  time  histories  in  Figure  9.  To  determine  these  equilibrium 
values  we  have  analyzed  the  cumulative  distributions  of  the 
Hj“Q— N2—O3  (*  =  5~"7)  torsional  angles  for  all  the  molecules 
in  the  simulation  box  during  trajectories  executed  at  different 
pressures,  T  =  293  K.  The  corresponding  data  obtained  in  these 
analyses  are  represented  in  Figure  4b.  When  compared  to  the 
distributions  of  these  angles  as  a  function  of  temperature  (see 
Figure  4a),  the  increase  of  pressure  has  a  totally  different  effect 
on  the  distribution  of  torsional  angles.  In  this  case  there  is  a 
continuous  shift  of  the  peak  positions  with  pressure  such  that 
between  0.3  and  5.4  GPa  this  shift  amounts  to  about  41°  while 
between  0.3  GPa  and  7.0  GP  the  corresponding  variation  is 
about  50°.  These  data  agree  very  well  with  the  experimental 
findings  obtained  by  Cromer  et  al.9  in  which  the  methyl  group 
was  found  to  be  rotated  by  about  45°  relative  to  die  low- 
temperature  configuration. 

We  next  explore  the  basis  for  the  rotational  reorientation  with 
pressure  by  comparing  properties  of  two  crystals,  one  in  which 
the  methyl  groups  are  rotated  and  one  in  which  the  methyl 
groups  are  in  the  arrangement  corresponding  to  the  low- 


temperature,  low-pressure  configuration.  Both  crystals  have 
P2i2i2i  symmetry  and  use  the  structural  parameters  correspond¬ 
ing  to  averages  obtained  in  our  simulations.  The  only  difference 
between  the  two  models  is  in  the  arrangement  of  the  hydrogens 
of  the  methyl  group.  In  one  model,  the  methyl  group  is  rotated 
by  50°  relative  to  the  low-temperature,  low-pressure  crystal.  In 
the  second  model,  the  methyl  group  is  in  the  same  arrangement 
as  the  low-temperature,  low-pressure  crystal.  We  first  compared 
the  total  potential  energy  for  the  two  systems  using  structural 
parameters  for  T  =  4.2  K  at  0  atm  pressure.  The  model  in  which 
the  methyl  group  is  not  rotated  has  a  total  potential  energy  that 
is  lower  by  3.5  kJ/mol  per  molecule.  In  that  model,  the  repulsive 
N— H  contributions  to  the  total  potential  are  greater  (by  4.9  kJ/ 
mol),  but  the  repulsive  H— H  interactions  are  lower  by  1.4  kJ / 
mol.  The  attractive  C-H  contribution  is  larger  for  the  model 
in  which  the  methyl  group  is  rotated  by  0.9  kJ/mol.  The  largest 
difference  in  contribution  to  the  total  potential  is  due  to  the 
O— H  contributions,  which  are  more  attractive  by  6.1  kJ/mol 
for  the  low-energy  crystal.  These  combined  effects  account  for 
the  difference  in  energy  between  the  two  models. 

For  model  crystals  using  structural  parameters  corresponding 
to  a  pressure  of  7  GPa  and  temperature  of  293  K,  the  model  in 
which  the  methyl  groups  are  rotated  by  50°  has  a  total  potential 
energy  that  is  lower  by  35  kJ/mol  per  molecule  than  the  model 
in  which  the  methyl  groups  are  not  rotated.  The  energy 
difference  (per  molecule)  is  due  mainly  to  the  H— N,  O— H,  and 
H-H  interactions.  For  the  model  in  which  the  methyl  group  is 
rotated  by  50°,  there  is  a  net  decrease  of  17.0  kJ/mol  for  the 
N— H  repulsive  interactions,  a  net  increase  of  1.1  kJ/mol  for 
the  C-H  attractive  interactions,  and  a  net  increase  of  26.8  kJ / 
mol  for  the  O-H  attractive  interactions.  H— H  repulsions  are 
greater  in  the  rotated  crystal  model  by  10  kJ/mol.  For  both 
pressures,  the  largest  differences  in  the  total  potential  energy 
are  due  to  the  contributions  of  O— H  and  N— H  interactions. 
The  results  indicate  that  the  methyl  group  rotation  with  pressure 
allows  for  the  enhancement  of  the  O—H  attractions  while 
reducing  the  N-H  repulsions. 

Structural  molecular  parameters  predicted  by  the  model  and 
experimental  values  (where  available)  as  a  function  of  pressure 
are  given  Table  9.  The  most  significant  changes  in  molecular 
structure  with  pressure  are  the  compression  of  the  C— N  bond 
and  the  distortion  of  the  Q—  N2—O4  angle.  These  trends  were 
also  seen  in  experiment;  however,  the  effects  were  more 
pronounced.  In  the  experimental  results,  the  C— N  bond  is 
compressed  by  0.08  A  at  3.5  GPa,  and  the  heavy  atom  angles 
have  changed  by  ~4°.  In  our  model,  the  corresponding  C— N 
bond  compression  is  only  0.01  A,  and  the  heavy  atom  angles 
have  not  significantly  changed. 

A  comparison  of  the  average  fractional  coordinates  and 
orientational  Euler  parameters  for  the  four  molecules  in  the  unit 
ceil  is  given  in  Table  9.  All  orientational  parameters  of  the  four 
molecules  in  the  unit  cell  were  averaged  over  time  and  all  unit 
cells  in  the  simulation  space.  We  have  also  generated  orienta¬ 
tional  parameters  for  the  four  molecules  in  the  unit  cell  assuming 
perfect  P2i2i2i  symmetry  in  the  manner  described  earlier.  The 
only  complete  set  of  fractional  coordinates  for  all  atoms  in  the 
Cromer  et  al.9  study  corresponded  to  3.5  GPa,  293  K,  with  the 
assumption  that  the  C— H  bonds  and  H— H  distances  are  1.0 
and  1.63  A,  respectively.  The  simulation  results  are  in  good 
agreement  with  these  experimental  data,  indicating  that  the 
model  reasonably  represents  the  crystal  under  moderate  pressure. 
Comparison  of  the  averages  of  the  four  molecules  in  the  unit 
cell  with  those  with  P2i2i2i  symmetry  indicate  that  the  space 
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TABLE  9:  Comparison  of  the  Average  Orientational  Parameters  (fractional  coordinates  sx,  sy,  sz  and  Euler  angles  0,  O,  V) 
for  the  Four  Molecules  in  the  Unit  Cell  and  of  the  Internal  Geometrical  Parameters  of  Nitromethane  Molecule  as  Obtained 
from  Trajectory  Calculations  at  Different  Pressures;  Corresponding  Experimental  Values  (where  available)  Are  Also  Indicated* 


parameter 

pressure(GPa) 

0.3 

NPT-MD 

0.3 

ideal 

3.5 

expt*^ 

3.5 

NPT-MD 

3.5 

ideal 

7.0 

NPT-MD 

7.0 

ideal 

sxl 

0.3488 

0.3453 

0.3391 

0.3316 

syl 

0.4284 

0.4055 

0.4311 

0.4304 

szl 

0.3765 

0.3844 

0.3855 

0.3938 

sx2 

0.1519 

0.1512 

0.1547 

0.1609 

0.1609 

0.1684 

0.1684 

sy2 

-0.4243 

-0.4284 

-0.4055 

-0.4300 

-0.4311 

-0.4295 

-0.4304 

sz2 

-0.1235 

-0.1235 

-0.1156 

-0.1139 

-0.1145 

-0.1059 

-0.1062 

sx3 

-0.3479 

-0.3488 

-0.3453 

-0.3388 

-0.3391 

-0.3314 

-0.3316 

sy3 

-0.0687 

-0.0716 

-0.0945 

-0.0688 

-0.0689 

-0.0693 

-0.0696 

sz3 

0.1233 

0.1235 

0.1156 

0.1149 

0.1145 

0.1060 

0.1062 

sx4 

-0.1525 

-0.1512 

-0.1547 

-0.1614 

-0.1609 

-0.1685 

-0.1684 

sy4 

0.0686 

0.0716 

0.0945 

0.0700 

0.0689 

0.0706 

0.0696 

sz4 

-0.3769 

-0.3765 

-0.3844 

-0.3853 

-0.3855 

-0.3940 

-0.3938 

01 

91.4 

91.3 

90.6 

89.5 

0)1 

54.7 

55.8 

56.2 

57.0 

'FI 

245.3 

242.9 

248.0 

249.3 

©2 

-91.4 

-91.4 

-91.3 

-89.9 

-90.6 

-89.1 

-89.5 

$2 

54.7 

54.7 

55.8 

56.3 

56.2 

57.0 

57.0 

W2 

245.3 

245.3 

242.9 

247.7 

248.0 

249.0 

249.3 

03 

88.4 

88.6 

88.7 

89.6 

89.4 

90.7 

90.5 

$3 

-54.8 

-54.7 

-55.8 

-56.3 

-56.2 

-57.1 

-57.0 

W3 

65.2 

65.3 

62.9 

67.8 

68.0 

69.0 

69.3 

©4 

-88.6 

-88.6 

-88.7 

-89.9 

-89.4 

-90.7 

-90.5 

<I>4 

-54.6 

-54.7 

-55.8 

-56.2 

-56.2 

-57.0 

-57.0 

65.2 

65.3 

62.9 

67.7 

68.0 

69.0 

69.3 

pressure(GPa) 

0.3 

0.3 

3.5 

3.5 

7.0 

parameter 

expt* 

NPT-MD 

expt* 

NPT-MD 

NPT-MD 

C-N 

1.55 

1.500 

1.47 

1.489 

1.481 

C-H 

1.092 

1.088 

1.084 

n-o3 

1.21 

1.231 

1.21 

1.228 

1.224 

n-o4 

1.20 

1.231 

1.26 

1.228 

1.226 

HCH 

110.5 

110.7 

110.6 

110.8 

HCN 

108.1 

107.9 

107.9 

107.7 

cno3 

117.7 

117.7 

118.4 

117.0 

115.9 

cno4 

116.2 

117.9 

119.9 

118.8 

119.9 

C>3N04 

125.9 

124.1 

.  121.6 

123.9 

123.9 

a  In  the  case  of  molecules  2—4,  the  ideal  orientational  parameters  obtained  based  on  P2i2j2i  space  group  symmetry  operators  and  the  orientation 
parameters  of  molecule  1,  positioned  in  the  asymmetric  unit,  are  also  indicated.  b  Experimental  values  from  ref  9.  Since  the  hydrogen  atom  positions 
could  not  be  resolved  at  0.3  GPa,  the  orientational  parameters  for  nitromethane  at  this  pressure  are  not  available  for  comparison.  c  Unit  cell  fractional 
values  have  been  shifted  to  fall  within  the  range  from  —1/2  to  1/2. 


group  symmetry  is  conserved  with  pressure,  as  observed  in 
experiment.9 

Also,  as  opposed  to  the  temperature  effects,  which  tend  to 
increase  the  degree  of  disorder  in  system,  the  increase  of 
pressure  has  the  opposite  effect.  As  can  be  seen  in  Figure  3c 
the  peaks  of  RDFs  for  O-C  intermolecular  distances  shift 
toward  smaller  distance  values  with  increasing  pressure,  indicat¬ 
ing  the  compression  of  the  material.  Moreover,  more  and  more 
peaks  regain  their  individuality  with  pressure  increases,  indicat¬ 
ing  an  increase  in  the  degree  of  rotational  and  translational  order. 

V.  Conclusions 

We  have  developed  a  classical  potential  for  simulation  of 
solid  nitromethane  containing  both  intra-  and  intermolecular 
potential  terms.  The  parametrization  of  the  intramolecular 
potential  has  been  done  on  the  basis  of  the  results  of  ab  initio 
calculations  performed  on  the  isolated  nitromethane  molecule 
at  the  B3LYP/6-31G*  level.  Both  structural  geometrical  pa¬ 
rameters  as  well  as  data  about  the  vibrational  frequencies  and 
their  eigenvectors  have  been  used  in  the  fitting  procedure.  The 


intermolecular  potential  used  in  these  calculations  is  of  the 
Buckingham  6-exp  form  and  was  previously  developed  for  RDX 
crystals1  and  showed  to  be  transferable  to  30  nitramine  crystals4 
and  to  5 1  other  non-nitramine  crystals.5 

Molecular  packing  calculations  using  the  proposed  set  of 
intermolecular  parameters  and  the  set  of  HF  charges  indicate 
an  accurate  prediction  of  crystallographic  parameters  with 
deviations  less  than  1.21%  for  the  lattice  edges.  Additionally, 
the  predicted  lattice  energy  of  nitromethane  is  in  acceptable 
agreement  with  previous  theoretical  estimations.34 

The  tests  of  this  potential  in  NPT-MD  simulations  indicate 
that  the  prediction  of  the  crystallographic  parameters  is  also 
well  reproduced  as  a  function  of  temperature  with  deviations 
between  1.7  and  3.6%  from  the  experimental  data7’9  for  the 
temperature  range  4.2—228  K.  Moreover,  the  present  potential 
is  able  to  predict  changes  of  the  crystallographic  parameters 
with  pressure  similar  to  those  observed  experimentally.9  A 
similar  good  agreement  is  found  for  the  bulk  modulus  and  its 
pressure  derivative.  At  zero  pressure  the  predicted  bulk  modulus 
is  Bo  =  6.78  GPa  while  the  corresponding  experimental  value 
is  7.0  GPa.9 


Solid  Nitromethane 

The  results  of  MD  and  NPT-MD  simulations  performed  using 
the  rigid-body  approximation  indicate  that  the  predicted  values 
in  these  cases  are  similar  to  those  obtained  from  the  flexible 
model  particularly  for  pressures  up  to  about  5  GPa. 

Both  RDFs  as  well  as  the  analyses  of  the  HCNO  dihedral 
angle  distributions  during  trajectories  calculations  obtained  using 
the  flexible  model  indicate  two  distinctive  regimes  for  methyl 
group  rotational  dynamics.  There  is  an  increase  of  the  degree 
of  rotation  of  the  methyl  group  with  temperature  increase  (at 
ambient  pressure)  without  significant  changes  of  the  average 
equilibrium  positions.  On  the  other  hand,  at  ambient  temperature, 
the  increase  of  pressure  causes  a  continuous  shift  of  the  average 
equilibrium  positions  of  the  H  atoms  relative  to  the  C-NO2 
plane.  A  comparison  of  the  contributions  of  the  individual 
intermolecular  interactions  to  the  total  potential  energies  for 
crystals  in  which  the  methyl  groups  are  in  the  low-temperature, 
low-pressure  configuration  or  arranged  in  the  high-pressure 
configuration  show  that  at  high  pressures,  the  O-H  attractive 
interactions  are  enhanced  upon  methyl  group  rotation  in  the 
crystal  while  the  N— H  repulsive  interactions  are  decreased;  the 
remaining  intermolecular  interactions  do  not  change  significandy 
between  the  two  models. 

The  success  of  the  present  potential  model  to  describe  the 
prototypical  explosive,  nitromethane,  conjugated  with  the 
performances  of  the  present  intermolecular  potential  to  describe 
a  large  number  of  important  energetic  materials,  provides  further 
incentives  to  develop  these  models  and  to  extend  their  applica¬ 
tions  to  more  subtle  effects  such  as  the  energy  transfer3738  and 
reactions  in  condensed  phases. 

Acknowledgment  We  are  pleased  to  acknowledge  many 
inspiring  and  helpful  discussions  with  Dr.  S.  F.  Trevino.  This 
work  was  supported  by  the  Strategic  Environmental  Research 
and  Development  Program  (SERDP).  D.L.T.  gratefully  ac¬ 
knowledges  support  by  the  U.S.  Army  Research  Office  under 
Grant  DAAG55-98- 1-0089. 

References  and  Notes 

(1)  Sorescu,  D.  C.;  Rice,  B.  M.;  Thompson,  D.  L.  J.  Phys.  Chem.  1997, 
B101,  798. 

(2)  Sorescu,  D.  C;  Rice,  B.  M.;  Thompson,  D.  L.  J.  Phys.  Chem.  1998, 
B102 ,  948. 

(3)  Sorescu,  D.  C.;  Rice,  B.  M.;  Thompson,  D.  L.  J.  Phys.  Chem.  1998, 
B  102 , 6692. 

(4)  Sorescu,  D.  C.;  Rice,  B.  M.;  Thompson,  D.  L.  J.  Phys.  Chem.  1998, 
A  102,  8386. 

(5)  Sorescu,  D.  C;  Rice,  B.  M.;  Thompson,  D.  L.  J.  Phys.  Chem.  1999, 
A  103 , 989. 

(6)  Sorescu,  D.  C.;  Rice,  B.  M.;  Thompson,  D.  L.  J.  Phys.  Chem.  1999, 
B  103 , 6783. 

(7)  Trevino,  S.  F.;  Prince,  E.;  Hubbard,  C.  R.  J.  Chem.  Phys.  1980, 
73 , 2996. 

(8)  Trevino,  S.  F.;  Rymes,  W.  H.  J.  Chem.  Phys.  1980,  75,  3001. 

(9)  Cromer,  D.  T.;  Ryan,  R.  R.;  Schiferl,  D.  J.  Phys.  Chem.  1985,  89, 
2315. 

(10)  GroSev,  V.  M.;  Stelzer,  F.;  Jocham,  D.  J.  Mol  Structure  1999, 476, 
181. 


J.  Phys.  Chem.  B,  Vol.  104 ,  No.  35,  2000  8419 

(11)  Jones,  W.  M.;  Giauque,  W.  F.  J.  Am.  Chem.  Soc.  1947,  69,  983. 

(12)  Yarger,  F.  L.;  Olinger,  B.  J.  Chem.  Phys.  1986,  85,  1534. 

(13)  Piermarini,  G.  J.;  Block,  S.;  Miller,  P.  J.  Phys.  Chem.  1989,  93, 
457. 

(14)  Cavagnat,  D.;  Mageri,  A.;  Vettier,  C.;  Anderson,  I.  S.;  Trevino,  S. 
F.  Phys.  Rev.  Lett.  1985,  54,  193. 

(15)  Tannebaum,  E.;  Myers,  R.  J;  Gwinn,  W.  D.  J.  Chem.  Phys.  1956, 
25,  42. 

(16)  Remizov,  A.  B.;  Musayakaeva,  R.  H.  Opt.  Spektrosk.  1975,  38, 
226. 

(17)  Frisch,  M.  J.;  Trucks,  G.  W.;  Schlegel,  H.  B.;  Gill,  P.  M.  W. 
Johnson,  B.  G.;  Robb,  M.  A.;  Cheeseman,  J.  R.;  Keith,  T.;  Patersson,  G. 

A. ;  Montgomery,  J.  A.;  Raghavachari,  K.;  Al-Laham,  M.  A.;  Zakrzewski, 

V.  G.;  Ortiz,  J.  V.;  Foresman,  J.  B.;  Cioslowski,  J.;  Stefanov,  B.  B.; 
Nanyakkara,  A.;  Challacombe,  M.;  Peng,  C.  Y.;  Ayala,  P.  Y.;  Chen,  W.; 
Wong,  M.  W.;  Andres,  J.  L.;  Replogle,  E.  S.;  Gomperts,  R.;  Martin,  R.  L.; 
Fox,  D.  J.;  Binkley,  J.  S.;  Defrees,  D.  J.;  Baker,  J.;  Stewart,  J.  P.;  Head- 
Gordon,  ML;  Gonzales,  C.;  Pople,  J.  A.  Gaussian  94,  Revision  C.3;  Gaussian, 
Inc.:  Pittsburgh,  PA,  1995. 

(18)  MoUer,  C.  M.  S.  Phys.  Rev .  1934,  46,  618. 

(19)  Rice,  B.  M.;  Thompson,  D.  L.  J.  Chem.  Phys.  1990,  93,  7986. 

(20)  (a)  Becke,  A.  D.  J.  Chem.  Phys.  1993, 98, 5648.  (b)  Lee,  C.;  Yang, 

W. ;  Pair,  R.  G.  Phys.  Rev.  1988,  B41,  785. 

(21)  Hehre,  W.;  Radom,  L.;  Schleyer,  P.  v.  R.;  Pople,  J.  A.  Ab  Initio 
Molecular  Orbital  Theory ;  Wiley-Interscience:  New  York,  1986. 

(22)  Rice,  B.  M.;  Chabalowski,  C.  F.;  Adams,  G.  F.;  Mowrey,  R.  C.; 
Page,  M.  Chem.  Phys.  Lett.  1991,  184,  335. 

(23)  Forseman,  J.  B.;  Frisch,  A.  In  Exploring  Chemistry  with  Electronic 
Structure  Methods',  Gaussian,  Inc.:  Pittsburgh,  PA,  1996. 

(24)  (a)  Wells,  A.  J.;  Wilson,  E.  B.  J.  Chem.  Phys.  1941,  9,  314.  (b) 
Smith,  D.  C.;  Pan,  C.;  Nielsen,  J.  R.  J.  Chem.  Phys.  1950,  18,  706.  (c) 
Jones,  W.  J.;  Sheppard,  N.  Proc.  R.  Soc.  London  Ser.  A  1968,  304,  135. 
(d)  Trinquecoste,  C.;  Rey-Lafon,  M.;  Forel,  M.-T.  Spectrochim.  Acta  A 
1974  30,  813.  (e)  McKean,  D.  C.;  Watt,  R.  A.  J.  Mol.  Spectrosc.  1976,  61, 
184. 

(25)  Pertsin,  A,  J.;  Kitaigorodsky,  A.  L  The  Atom-Atom  Potential 
Method,  Applications  to  Organic  Molecular  Solids',  Springer-Verlag:  Berlin, 
1987. 

(26)  Desiraju,  G.  R.  Crystal  Engineering:  The  Design  of  Organic  Solids; 
Elsevier  Amsterdam,  1989. 

(27)  Gibson,  K.  D.;  Scheraga,  H.  A.  J.  Phys.  Chem.  1995,  99,  3752. 

(28)  Gibson,  K.  D.;  Scheraga,  H.  A.  LMIN:  A  Program  for  Crystal 
Packing,  QCPE,  No.  664. 

(29)  Melchionna,  S.;  Ciccotti,  G.;  Holian,  B.  L.  Mol.  Physics  1993,  78, 
533. 

(30)  DL_POLY  is  a  package  of  molecular  simulation  routines  written 
by  W.  Smith  and  T.  R.  Forester,  copyright  The  Council  for  the  Central 
Laboratory  of  the  Research  Councils,  Daresbury  Laboratory  at  Daresbury, 
Nr.  Warrington,  1996. 

(31)  Allen,  M.  P.;  Tildesley,  D.  J.  Computer  Simulation  of  Liquids; 
Oxford  University  Press:  New  York,  1989. 

(32)  Fincham,  D.  Molecular  Simulation  1992,  8,  165. 

(33)  Fundamentals  of  Crystallography;  Giacovazzo,  C.,  Ed.;  Oxford 
University  Press:  New  York,  1992. 

(34)  The  enthalpy  of  sublimation  for  nitromethane  has  been  com¬ 
municated  by  P.  Politzer  and  has  been  obtained  based  on  the  method 
described  in  Politzer,  P.;  Murray,  J.  S.;  Grice,  M.  E.;  DeSalvo,  M.;  Miller, 
E.  M.  Mol.  Phys.  1997,  92,  923. 

(35)  Mumaghan,  F.  D.  In  Finite  Deformation  of  an  Elastic  Solid;  Dover 
Publications:  New  York,  1951;  p  73. 

(36)  Cook,  M.  D.;  Fellows,  J.;  Haskins,  P.  J.  In  Decomposition, 
Combustion  and  Detonation  Chemistry  of  Energetic  Materials;  Brill,  T. 

B. ,  Russell,  T.  P.,  Tao,  W.  C.,  Wardle,  R.  B.,  Eds.;  Mater.  Res.  Soc.  Symp. 
Proc .  1995,  418,  267-275. 

(37)  Aubuchon,  C.  M.;  Rector,  K.  D.;  Holmes,  W.;  Fayer,  M.  D.  Chem. 
Phys.  Lett.  1999,  200,  84. 

(38)  Deak,  J.  C.;  Iwaki,  L.  K.;  Dlott,  D.  D.  J.  Phys.  Chem.  A  1999, 
103,  971. 


Intentionally  left  blank. 


NO.  OF 
COPIES 

2 


1 


1 


1 


I 


1 


1 


NO.  OF 

ORGANIZATION  COPIES 

DEFENSE  TECHNICAL 

INFORMATION  CENTER  1 

DTIC  DDA 

8725  JOHN  J  KINGMAN  RD 
STE  0944 

FT  BELVOIR  VA  22060-6218 

HQDA  1 

DAMOFDT 

400  ARMY  PENTAGON 

WASHINGTON  DC  20310-0460 

OSD 

OUSD(A&T)/ODDDR&E(R)  3 

RJTREW 

THE  PENTAGON 

WASHINGTON  DC  20301-7100 

DPTY  CG  FOR  RDA 

US  ARMY  MATERIEL  CMD  1 

AMCRDA 

5001  EISENHOWER  AVE 
ALEXANDRIA  VA  22333-0001 

INST  FOR  ADVNCD  TCHNLGY 
THE  UNTV  OF  TEXAS  AT  AUSTIN 
PO  BOX  202797 
AUSTIN  TX  78720-2797 

4 

DARPA 

BKASPAR 

3701  N  FAIRFAX  DR 

ARLINGTON  VA  22203-1714 

US  MILITARY  ACADEMY 

MATH  SCI  CTR  OF  EXCELLENCE 

MADN  MATH 

MAJ  HUBER 

THAYER  HALL 

WEST  POINT  NY  10996-1786 

DIRECTOR 

US  ARMY  RESEARCH  LAB 

AMSRLD 

DR  SMITH 

2800  POWDER  MILL  RD 
ADELPHI MD  20783-1 197 


ORGANIZATION 


DIRECTOR 

US  ARMY  RESEARCH  LAB 
AMSRLDD 

2800  POWDER  MILL  RD 
ADELPHI  MD  20783-1 197 

DIRECTOR 

US  ARMY  RESEARCH  LAB 
AMSRL  Q  AIR  (RECORDS  MGMT) 
2800  POWDER  MILL  RD 
ADELPHI  MD  20783-1 145 

DIRECTOR 

US  ARMY  RESEARCH  LAB 
AMSRL  Cl  LL 
2800  POWDER  MILL  RD 
ADELPHI  MD  20783-1145 

DIRECTOR 

US  ARMY  RESEARCH  LAB 
AMSRL  a  AP 
2800  POWDER  MOLL  RD 
ADELPHI  MD  20783-1 197 


ABERDEEN  PROVING  GROUND 
DIRUSARL 

AMSRL  a  LP  (BLDG  305) 


NO.  OF 
COPIES 


ORGANIZATION 


ABERDEEN  PROVING  GROUND 

22  DIRUSARL 
AMSRLWM 
B  RINGERS 
AMSRLWM  BD 
BE  PORCH 
WR  ANDERSON 
SWBUNTE 
C  F  CHABALOWSK3 
A COHEN 
R  DANIEL 
D  DEVYNCK 
R  AFDFER 
BE  HOMAN 
AJKOTLAR 
K  L  MCNESBY 
M  MCQUAID 
MS  MILLER 
A  W  MIZIOLEK 
JB  MORRIS 

R  A  PESCE-RODRIGUEZ 
BMRICE 
R  C  SAUSA 
MASCHROEDER 
J  A  VANDERHOFF 
AMSRLWM  MB 
B  FINK 


REPORT  DOCUMENTATION  PAGE 

Form  Approved 

OMB  No.  0704-0188 

1 .  AGENCY  USE  ONLY  (iMIV*  blank)  2.  REPORT  DATE  3.  REPORT  TYPE  AND  DATES  COVERED 

February  2001  Reprint,  Jun  1999-Feb  2000 

4.  TITLE  AND  SUBTITLE 

Theoretical  Studies  of  Solid  Nitromethane 

5.  FUNDING  NUMBERS 

622618.H80 

6.  AUTHOR(S) 

Dan  C.  Sorescue,  *  Betsy  M.  Rice,  and  Donald  L.  Thompson* 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

U.S.  Army  Research  Laboratory 

ATTN:  AMSRL-WM-BD 

Aberdeen  Proving  Ground,  MD  21005-5066 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

ARL-RP-14 

9.  SPONSORING/M ONITORING  AGENCY  NAMES(S)  AND  ADDRESS(ES) 

1 0.SPONSORING/MONITORING 

AGENCY  REPORT  NUMBER 

11.  SUPPLEMENTARY  NOTES 

*  Department  of  Chemistry,  Oklahoma  State  University,  Stillwater,  OK  74078 

A  reprint  from  the  Journal  of  Physical  Chemistry  By  vol.  104,  no.  35,  pp.  8406-8419. 

12*.  DtSTRIBUTION/AVAi  LABILITY  STATEMENT 

Approved  for  public  release;  distribution  is  unlimited. 

12b.  DISTRIBUTION  CODE 

13.  ABSTRACT  (Maximum 200  wonty 

A  classical  potential  to  simulate  the  dynamics  of  a  nitromethane  crystal  as  a  function  of  temperature  and  pressure  is 
described.  The  intramolecular  part  of  the  potential  was  taken  as  superposition  of  bond  stretching,  bond  bending,  and 
torsional  angles  terms.  These  terms  were  parametrized  on  the  basis  of  the  geometric  and  spectroscopic  (vibrational 
frequencies  and  eigenvectors)  data  obtained  using  ab  initio  molecular  orbital  calculations  performed  at  the 
B3LYP/6-31G  level  on  an  isolated  molecule.  The  intermolecular  potential  used  is  of  the  Buckingham  6-exp  form  plus 
charge-charge  Coulombic  interactions  and  has  been  previously  developed  by  us  (Sorescu,  D.  C.;  Rice,  B.  M.; 
Thompson,  D.  L.  J.  Phys.  Chem.  1997,  B101,  798)  to  simulate  crystals  containing  nitramine  molecules  and  several 
other  classes  of  nitro  compounds.  The  analyses  performed  using  constant  pressure  and  temperature  molecular 
dynamics  simulations  and  molecular  packing  calculations  indicate  that  the  proposed  potential  model  is  able  to 
reproduce  accurately  the  changes  of  the  structural  crystallographic  parameters  as  functions  of  temperature  or  pressure 
for  the  entire  range  of  values  investigated.  In  addition,  the  calculated  bulk  modulus  of  nitromethane  was  found  in 
excellent  agreement  with  the  corresponding  experimental  results.  Moreover,  it  was  determined  that  the  present 
potential  predicts  correctly  an  experimentally  observed  45°  change  in  methyl  group  orientation  in  the  high-pressure 
regime  relative  to  the  low-temperature  configuration. 

14.  SUBJECT  TERMS 

nitromethane,  molecular  dynamics,  intermolecular  interaction  potential 

15.  NUMBER  OF  PAGES 

20 

16.  PRICE  CODE 

17.  SECURITY  CLASSIFICATION  18.  SECURITY  CLASSIFICATION  19.  SECURITY  CLASSIFICATION 

OF  REPORT  OF  THIS  PAGE  OF  ABSTRACT 

UNCLASSIFIED  UNCLASSIFIED  UNCLASSIFIED 

20.  LIMITATION  OF  ABSTRACT 

UL 

NSN  7540-01-280-5500  Standard  Form  298  (Rev.  2-89) 


Prescribed  by  ANSI  Std.  239-18  298-102 


Intentionally  left  blank. 


USER  EVALUATION  SHEET/CHANGE  OF  ADDRESS 


This  Laboratory  undertakes  a  continuing  effort  to  improve  the  quality  of  the  reports  it  publishes.  Your  comments/answers  to 
the  items/questions  below  will  aid  us  in  our  efforts. 

1.  ARL  Report  Number/Author  ARL-RP-14  (POC:  Rice) _ Date  of  Report  February  2001 _ 

2.  Date  Report  Received _ _ _ _ _ _ _ _ 

3.  Does  this  report  satisfy  a  need?  (Comment  on  purpose,  related  project,  or  other  area  of  interest  for  which  the  report  will  be 

used.)  _ _ _ _ _ _ 


4.  Specifically,  how  is  the  report  being  used?  (Information  source,  design  data,  procedure,  source  of  ideas,  etc.) 


5.  Has  the  information  in  this  report  led  to  any  quantitative  savings  as  far  as  man-hours  or  dollars  saved,  operating  costs 
avoided,  or  efficiencies  achieved,  etc?  If  so,  please  elaborate.__ _ 


6.  General  Comments.  What  do  you  think  should  be  changed  to  improve  future  reports?  (Indicate  changes  to  organization, 
technical  content,  format,  etc.) _ _ _ _ _ 


Organization 

CURRENT  Name  E-mail  Name 

ADDRESS  _ 

Street  or  P.O.  Box  No. 

City,  State,  Zip  Code 

7.  If  indicating  a  Change  of  Address  or  Address  Correction,  please  provide  the  Current  or  Correct  address  above  and  the  Old  or 
Incorrect  address  below. 


Organization 


OLD  Name 

ADDRESS  _ 

Street  or  P.O.  Box  No. 


City,  State,  Zip  Code 
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